SCATTERING  MATRIX  ELEMENTS  FOR  THE 
NONADIABATIC  COLLISION 


DISSERTATION 


Luke  A.  Barger,  Captain,  USAF 
AFIT/DS/ENP/ 1 0-S02 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

Air  Force  Institute  of  Technology 

Wright-Patterson  Air  Force  Base,  Ohio 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

01  DEC  2010 


2.  REPORT  TYPE 


3.  DATES  COVERED 

00-00-2010  to  00-00-2010 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


4.  TITLE  AND  SUBTITLE 

SCATTERING  MATRIX  ELEMENTS  FOR  THE 


6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING  ORGANIZATION 

Air  Force  Institute  of  Technology, 2950  Hobson  report  number 

Way, WPAFB, OH, 45433-7765 

9.  SPONSORING/MONITORING  AGENCY  NAME(S )  AND  ADDRESS(ES )  10.  SPONSOR/MONITOR' S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

Scattering  matrix  elements  are  calculated  for  the  nonadiabatic  inelastic  collision  .  This  calculation  utilizes 
the  effective  potential  energy  surfaces  for  this  collision  generated  by  Garvin  along  with  a  correction  to  the 
asymptotic  H2  potential.  Wavepackets  are  propagated  on  these  surfaces  using  a  split-operator  propagator. 
This  propagation  yields  correlation  functions  between  reactant  and  product  M?ller  states  which  are  used 
to  calculate  the  scattering  matrix  elements  with  the  channel  packet  method.  These  scattering  matrix 
elements  represent  probability  amplitudes  for  the  collision  to  result  in  changes  to  the  electronic  fine 
structure  and  to  the  rotational  and  vibrational  eigenstates  of  the  H2  molecule  over  a  range  of  energies,  and 
are  presented,  discussed  and  compared  to  previous  work  in  which  the  hydrogen  bond  length  was  fixed  at 
its  equilibrium  value.  A  method  for  approximating  probability  for  the  reaction  B+H2-BH+H  as  a  function 
of  collisional  energy  is  presented. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

132 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


The  views  expressed  in  this  dissertation  are  those  of  the  author  and  do  not  reflect  the 
official  policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense,  or  the 
United  States  Government. 


AFIT/DS/ENP/ 1 0-S02 


SCATTERING  MATRIX  ELEMENTS  FOR  THE  NONADIABATIC  COLLISION 


Bilpj,)+MX 


DISSERTATION 


Presented  to  the  Faculty 

Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
In  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Doctor  of  Philosophy 


Luke  A.  Barger,  BSE,  MS 
Captain,  USAF 

December  2010 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AFIT/DS/ENP/ 1 0-S02 


SCATTERING  MATRIX  ELEMENTS  FOR  THE  NONADIABATIC  COLLISION 


Luke  A.  Barger,  BSE,  MS 
Captain,  USAF 


Approved: 


_ Signed _  _ 

David  E.  Weeks,  PhD  (Chairman)  Date 


Signed _  _ 

Kevin  C.  Gross,  PhD  (Member)  Date 


_ Signed _  _ 

Michael  J.  Havrilla,  PhD  (Member)  Date 


Accepted: 


_ Signed _  _ 

M.  U.  Thomas  Date 

Dean,  Graduate  School  of 
Engineering  and  Management 


AFIT/DS/ENP/ 1 0-S02 


Abstract 

Scattering  matrix  elements  are  calculated  for  the  nonadiabatic  inelastic  collision 

B  ( 2Pj  )  +  //2  ( 1 ,  v,  y )  i?  1 j  +  H2  ( 1 £* ,  v j ') .  This  calculation  utilizes  the 

effective  potential  energy  surfaces  for  this  collision  generated  by  Garvin  [1]  along  with  a 
correction  to  the  asymptotic  H2  potential.  Wavepackets  are  propagated  on  these  surfaces 

using  a  split-operator  propagator.  This  propagation  yields  correlation  functions  between 
reactant  and  product  Moller  states  which  are  used  to  calculate  the  scattering  matrix 
elements  with  the  channel  packet  method  [2],  These  scattering  matrix  elements  represent 
probability  amplitudes  for  the  collision  to  result  in  changes  to  the  electronic  fine  structure 
and  to  the  rotational  and  vibrational  eigenstates  of  the  H2  molecule  over  a  range  of 

energies,  and  are  presented  and  discussed.  They  are  also  compared  to  previous  work  in 
which  the  hydrogen  bond  length  was  fixed  at  its  equilibrium  value.  A  method  for 
approximating  probability  for  the  reaction  B  +  H2—>  BH  +  H  as  a  function  of  collisional 
energy  is  presented. 
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SCATTERING  MATRIX  ELEMENTS  FOR  THE  NONADIABATIC  COLLISION 


I.  Introduction 


The  B+H2  System 

The  B  +  H2  system  became  of  particular  interest  when  boron  was  proposed  as  an 

energetic  additive  to  solid  molecular  hydrogen  (SMH)  for  use  as  a  novel  rocket 
propellant.  Realization  of  this  technology  would  require  complete  characterization  of  the 
chemical,  physical,  and  transport  properties  of  boron  doped  SMH  [3].  These  early 
studies  prompted  additional  theoretical  interest  in  the  B  +  H2  system.  Specifically 
challenging,  and  hence  interesting,  is  the  unpaired  electron  in  the  2  p  orbital  of  atomic 

boron.  This  electron  creates  a  2P  spectroscopic  term  consisting  of  three  degenerate 
orbitals.  When  in  proximity  of  a  hydrogen  molecule  this  degeneracy  is  slightly  lifted 
resulting  in  three  orbitals  whose  energy  levels  are  very  similar.  As  a  result,  the  electronic 
eigenstates  exhibit  strong  dependence  of  nuclear  coordinates,  causing  the  system  to 
behave  nonadiabatically  [4].  Given  the  P  electronic  character  of  the  B  +  H2  system  [3] 
[5]  [6]  [7],  it  can  also  serve  as  an  accessible  means  to  understand  other  atomic  systems 
that  exhibit  2P  electronic  behavior.  Of  particular  interest  is  the  insight  the  B  +  H , 
system  lends  to  collisional  de-excitation  in  the  operation  of  diode-pumped  alkali  lasers 


1 


(DPALs).  The  Air  Force  is  currently  very  interested  in  the  development  of  DPAL 
systems. 

During  operation  of  a  DPAL  this  electron  is  pumped  to  the  2  PV2  states  are 
collisionally  de-excited  to  the  2  PU2  states.  This  process  is  shown  for  rubidium  [8] 


5  2p 

J  1  3/2 

5 2  p 

j  1 1/2 


5  5, 


1/2 


Figure  1.  DPAL  Transition  for  Rubidium 


For  practical  operation,  electron  population  excited  to  the  2  Pm  level  must 
transition  to  the 2  PV2  level  at  a  sufficient  rate.  Specifically  it  must  be  much  faster  than  the 

2P  spontaneous  emission  rate  of  :  3xl07  sec  1  [8].  Ground-state  boron  also  has  a  2P 
spectroscopic  term.  Since  it  shares  the  same  valence  electron  configuration  with  alkalis 
in  the  first  excited  state  it  will  behave  similarly.  The  fine-structure  collisional  de¬ 
excitation  from  2  PV2  to  2  PV2  is  among  the  transitions  examined  when  studying  the 
collision  of  boron  with  molecular  hydrogen  and  will  provide  insight  into  the  behavior  of 
alkali  atoms  in  DPAL  systems. 
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Statement  of  Objectives 

This  work  will  investigate  the  nonadiabatic  effects  of  the  collision  of  a  boron 
atom  with  molecular  hydrogen,  given  by  the  equation 

M  5pJ + "4 1  d  «  M  T;  J + "T  > ') 

Effective  potential  energy  surfaces  for  this  interaction  are  provided  by  Garvin  [1], 
extending  initial  ab  initio  calculations  performed  by  Dr.  David  Yarkony  at  Johns  Hopkins 
University. 

By  propagating  wavepackets  on  these  surfaces,  scattering  matrix  elements  can  be 
obtained  using  the  channel  packet  method  [2].  These  represent  energy-resolved 
probability  amplitudes  for  scattering  from  an  incoming  state  to  an  outgoing  state,  both  of 
which  are  selected  prior  to  the  calculation.  As  shown  by  eq.  (1),  the  collision  may 
involve  changes  to  the  hydrogen  vibrational  and/or  rotational  states  and  the  fine-structure 
of  the  boron  electronic  state.  The  objectives  of  this  research  are  to:  1)  characterize  the 
B  +  H2  effective  potential  energy  surfaces;  2)  propagate  wavepackets  on  these  surfaces; 

and  3)  compute  B  +  scattering  matrix  elements  and  present  analysis. 


(1) 
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II.  Theory 


Schrodinger’s  Equations 

Fundamental  to  much  of  quantum  mechanics,  and  this  research  in  particular,  is 
detennining  the  time-evolution  of  a  system.  The  governing  equation  is  the  time- 
dependent  Schrodinger  equation 

-nfkMk)  (2) 

where  |^)  is  the  wavefunction  and  H  is  the  Hamiltonian. 

The  time-dependent  Schrodinger  equation  is  most  easily  solved  when  the 
wavefunction  is  chosen  to  be  a  superposition  of  eigenstates  of  the  Hamiltonian 

i 

where  the  eigenstates  |  (f>)  are  solutions  to  the  equation 

«w=.Ek)  (4) 

Eq.  (4)  is  known  as  the  time-independent  Schrodinger  equation.  The  eigenvalue  E 
corresponding  to  each  eigenstate  is  the  energy  of  that  eigenstate.  The  set  of  all 
eigenstates  of  the  Hamiltonian  is  complete,  allowing  any  wavefunction  to  be  constructed 
as  a  superposition  of  eigenstates  of  the  Hamiltonian.  Solving  the  time-independent 
Schrodinger  equation  is  often  the  first  and  most  important  step  in  understanding  a 
quantum  mechanical  system. 
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The  Hamiltonian 


The  Hamiltonian  for  a  system  is  the  sum  of  the  kinetic  and  potential  energy 
operators  of  the  system, 

H  =T  +  V  (5) 

where  f  is  the  kinetic  energy  operator  and  V  is  the  potential  energy  operator.  If  the 
system  is  made  up  of  interacting  nuclei  and  electrons,  it  may  be  represented  as 

H  =  fN+f+Vm  +  VeN  +  Vee  (6) 

In  eq.  (6)  the  subscripts  identify  the  particle  types  -  N  for  nuclei  and  e  for  electrons.  The 
kinetic  energy  operators  are  simply  the  sum  of  the  kinetic  energy  operators  of  their 
respective  particle  type,  while  the  potential  operators  are  the  sum  of  the  Coulomb 
potentials  for  all  the  particles  -  among  nuclei,  between  nuclei  and  electrons,  and  among 
electrons,  respectively.  If  the  system  contains  nN  nuclei  and  ne  electrons,  the 
Hamiltonian  is 


A2  ”*  A2 

Pa  ,  V  Pi 


..,2 m„  2m,  „lj^qa-qf 


nN  ",  7  n2-  ",  ",  „?■ 

_yy^^+yy— — 

a= 1  != 1  Ha  ~  i=l  j>i  C1 ,  ~  j 


(7) 


where  p  is  the  momentum  operator,  in  is  the  mass  of  the  particle,  Ze  is  the  charge  of  the 
particle  (given  by  its  atomic  number  multiplied  by  the  charge  of  the  electron),  and  q  is 
the  position  operator  of  the  particle.  The  Greek  letters  a  and  P  index  nuclei,  while  i  and 
j  index  electrons. 

The  total  number  of  terms  in  the  Hamiltonian  increases  rapidly  as  the  number  of 
interacting  particles  increases.  The  B  +  H2  system  has  3  nuclei  and  7  electrons,  giving  a 


total  of  55  terms  in  the  Hamiltonian.  The  time-independent  Schrodinger  equation  for  the 


5 


Hamiltonian  given  in  eq.  (7)  cannot  be  solved  analytically,  prompting  the  use  of 
approximations  and  numerical  methods. 


The  Born-Oppenheimer  Approximation 

Since  nuclei  are  much  more  massive  than  electrons,  it  is  expected  that  they  will 
move  very  slowly  relative  to  the  electrons’  motion.  The  electrons  move  under  the 
influence  of  nuclei  which  are  approximately  fixed  in  space.  In  turn,  nuclei  move  under 
the  influence  of  an  average  electron  ‘cloud’.  This  suggests  an  approximation  for  solving 
the  molecular  time-independent  Schrodinger  equation  by  separating  nuclear  and 
electronic  motion. 

The  nuclei  are  fixed  in  place  by  setting  their  kinetic  energy  operators  to  zero  in 
the  Hamiltonian.  The  remaining  terms  are  known  as  the  electronic  Hamiltonian. 


EaE pe 


i= 1  2me 


-tl 


Ze 


-  + 


iii 


*e  a=\  p>a  yla  Q  p\  a=1  i=1  \^a  Qi  \  i= 1  j>i  \tfi  Q  j  | 

In  eq.  (8)  the  nuclear  coordinates  are  parameters  rather  than  operators,  and  the  second 


tenn  ^m(q,,)  =  XX 


"N  »J V  7  7  p 


is  a  constant  for  a  given  set  of  nuclear  coordinates 


a~\  P>a  <1  p\ 


N  {^1  ’  ^2  5  — }  ‘ 


The  electronic  time-independent  Schrodinger  equation  is 

Helec  |  j  {qN  ))  =  EJelec  (qN )  |  j(qN  )) 


(8) 


(9) 


where  both  the  electronic  eigenstates  |y)  and  their  associated  energy  eigenvalues  EJelec 
are  explicitly  shown  to  depend  parametrically  on  the  nuclear  coordinates  qN .  This 
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equation  may  be  solved  using  ab  initio  numerical  techniques  such  as  the  Hartree-Fock 
approximation,  configuration  interaction  or  multi-configuration  self-consistent  field 
methods.  These  techniques  are  described  by  Szabo  and  Ostlund  [9],  and  are  generally 
applied  through  the  use  of  mature  software  packages  developed  specifically  for  large 
quantum  chemistry  applications  to  yield  EJelec  over  a  range  of  nuclear  coordinates.  The 
result  is  a  potential  energy  surface  (PES)  EJelec  (qN) ,  which  has  a  dimensionality  given  by 
the  number  of  nuclear  degrees  of  freedom,  as  well  as  the  associated  electronic  eigenstates 
|./(9V))- 

The  set  of  electronic  eigenstates  is  complete  and  the  eigenstates  of  the  full  system 
|  y/ 'j  can  therefore  be  represented  in  the  basis  of  the  electronic  eigenstates. 

k(^))=ZF/(^)|y(^))  (io) 

j 

Substituting  this  into  the  full  time-independent  Schrodinger  equation  gives 

j  j 

This  expression  is  simplified  by  noting  that  E  is  a  constant  and  the  electron  eigenstates 
are  orthonormal.  Multiplying  from  the  left  by  (i (qN  ) |  and  integrating  over  electronic 
coordinates  gives 

{i(<l«)\{b+H^)Y,Ft  =  EF,(q„)  (12) 

j 

Additionally  the  electronic  Hamiltonian  does  not  operate  on  nuclear  coordinates, 
allowing  this  to  be  simplified  as 

{KVn  )|tl^/(9,)|./(^,v))  +  (*  fc  )  |  Z  Fi  fc  Pelec  \  J  {<1 N  )}  =  EF,  (<7,V  ) 

j  j 
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(13) 


K9»)|4Z^(9»)|;(9»))+^(9»)=^fc) 

j 

where  Helec\j{qN))  =  EJelec\j{qN))  and  ^i{qN)\j{qN))  =  5i].  The  nuclear  kinetic  energy 
operator  TN  does  operate  on  nuclear  coordinates.  Expressing  eq.  (13)  in  the  coordinate 

h2 

representation  with  |  and  TN  —>  ^ - -V2^  gives 

a  2tna 


(  )  VL  Fj  ( <7v  )  <t>j  (  <1n  >  Cle  )  F1  e  +  FLF ]  {<1 N)  =  FFi  {<1  n) 


The  dependence  of  the  various  functions  on  nuclear  coordinates  qN  and  electronic 
coordinates  qe  is  shown  explicitly.  The  subscript  V2  indicates  which  coordinates  the 
derivatives  are  with  respect  to.  Applying  the  product  rule  eq.  (14)  becomes 


(14) 


~  y  X  ~  K Vl  +  2R  ( 4  v ,  qe  ;  ( qN ,  qe )  dqe  ■  V?w 

Z  aj  ma 

+  q,  K «*,(«»,  9« ) dq.}  F,  ( q„  )  (15) 

=  (E~EL)F,(q„) 

This  expression  is  simplified  by  defining  two  terms.  The  first  is  a  vector  quantity  known 
as  the  derivative  coupling  tenn  (DCT): 

h  =  (16) 


The  second  is  a  scalar  quantity  known  as  the  kinetic  coupling  term: 

Kij  =  tJV*  {Qn’  Ve  )VL  ^  (  <iN  5  <le  )  dcle 


(17) 


Using  eqs.  (16)  and  (17)  gives 


TEl^-Kvl+2vi,+2*-d  os) 


2  77  ma 
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When  the  derivative  and  kinetic  coupling  terms  are  small  they  can  be  ignored,  and  the 
system  is  said  to  behave  adiabatically.  This  is  the  Bom-Oppenheimer  approximation. 
Eliminating  those  terms  from  eq.  (18)  gives  the  simpler  form 


I-' 

« 


Qn 


+  E 


elec 


FMn)  =  EFMn) 


(19) 


Here  eq.  (19)  has  been  cast  in  the  form  of  the  time-independent  Schrodinger  equation 
where  Ft  ( qN )  is  interpreted  as  the  nuclear  wavefunction.  This  wavefunction  only 

depends  on  one  adiabatic  PES,  E'elec .  If  the  derivative  and  kinetic  coupling  terms  cannot 

be  ignored,  the  system  is  said  to  behave  nonadiabatically,  and  the  coupling  in  eq.  (18) 
cannot  be  neglected.  For  this  reason  the  derivative  and  kinetic  coupling  terms  are 
collectively  known  as  nonadiabatic  coupling  terms  (NACTs). 


The  Adiabatic  to  Diabatic  Transformation 

Many  polyatomic  systems  behave  nonadiabatically  [10],  that  is,  they  contain 
nonnegligible  NACTs  when  expressed  in  the  adiabatic  representation  as  in  eq.  (18).  The 
clearest  predictor  of  nonadiabatic  behavior  is  when  the  adiabatic  electronic  PESs  E'elec  are 

close  together  [11].  This  is  seen  in  the  generalized  Hellman-Feynman  theorem  which 
gives  the  derivative  coupling  terms  in  the  adiabatic  representation  [12]: 


xj  = 


jM) 


FJ  —  Fl 

^ elec  ^ elec 


(20) 


The  B  +  H2  system  is  an  example  of  a  system  with  nonnegligible  NACTs.  The 


unpaired  2 p  electron  in  ground-state  boron  gives  rise  to  a  2  P  spectroscopic  tenn 
involving  three  degenerate  atomic  orbitals.  This  degeneracy  is  slightly  lifted  when  in  the 
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presence  of  a  hydrogen  molecule,  but  they  remain  energetically  very  close.  Eq.  (20) 
predicts  that  these  three  orbitals  will  be  strongly  coupled  to  each  other,  but  only  weakly 
coupled  to  the  other  atomic  orbitals  from  when  they  are  separated  energetically. 

The  problem  in  solving  the  coupled  set  of  equations  on  the  adiabatic  basis  is 
handled  by  transfonning  to  the  diabatic  basis  in  which  the  NACTs  are  negligible  [13]. 
This  is  accomplished  via  a  unitary  transfonnation  and  results  in  a  form  of  the  time- 
independent  Schrodinger  equation  where  the  wavefunction  is  coupled  by  the  diabatic 
PES. 

In  the  interest  of  computational  efficiency,  the  basis  of  electronic  eigenstates  is 
truncated  to  only  include  those  states  that  are  strongly  coupled.  These  are  the  three 
adiabatic  PESs  corresponding  to  the  three 2  P  states.  These  PESs  are  labeled  V2 ,  fk 

and  VA, ,  where  the  subscripts  are  inherited  from  the  Cs  symmetry  displayed  by  the 
B  +  H2  system  [1].  For  this  calculation  the  V2  ,  V2  and  VA„  PESs  were  calculated  by 


Dr.  David  Yarkony  using  a  multi-reference  configuration  interaction  calculation. 

The  adiabatic  time-independent  Schrodinger  equation  (18)  is  represented  in  the 
truncated  basis  consisting  of  these  three  states. 


0 

E 

0 


(21) 
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Since  the  basis  has  been  truncated  to  include  three  PES,  the  adiabatic  to  diabatic 
transformation  matrix  is  3x3.  The  restriction  that  the  electronic  wavefunctions  are  real, 


along  with  the  fact  that  states  with  A "  symmetry  type  will  not  mix  with  those  with  A ' 
symmetry  means  the  transfonnation  matrix  may  be  represented  by  an  orthogonal  matrix 


parameterized  by  a  single  angle  [14],  [15]: 


U ad  (*2v) 


"cos(y(^))  -sin(y(^)) 
sin(r(^))  cos(f(^)) 

0  0 


0^ 

0 


(22) 


This  transfonnation  matrix  relates  the  adiabatic  wavefunctions  to  the  diabatic 
wavefuctions  as 

Ka  (?» ><ie)  =  (?»  ’  Ve  )  sin  (r  {<ln  ))  +  4>*  A  (?„ »  <le  )  C0S  (/ (?„  )) 

Substituting  eq.  (23)  into  the  definition  of  the  derivative  coupling  terms  in  eq.  (16)  gives 
an  expression  for  the  derivative  coupling  terms  in  the  diabatic  representation  in  terms  of 
the  derivative  coupling  terms  in  the  adiabatic  representation  and  the  mixing  angle  y .  The 
derivative  coupling  terms  may  be  shown  to  be  antihermitian,  which  simplifies  this 
relation  to 

k=Kr(< i.)+k  m) 

The  requirement  that  the  mixing  angle  causes  the  diabatic  coupling  terms  to  be  zero  gives 

k  =-V,/(9„)  (25) 

and  the  mixing  angle  y(qN)  may  be  calculated  by  performing  a  line  integral  through  the 
vector  space  of  the  derivative  coupling  terms  at  each  nuclear  configuration: 
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(26) 


I\  (7  1 

r(R,0, r )  =  J  dR  +  J  r6»,i2 1 Rr^  dO  +  J rrl2 dr  +  y0 
*0  #0  r0 


The  adiabatic  to  diabatic  transfomiation  is  applied  to  the  truncated  basis 
representation  in  eq.  (21)  to  give  the  diabatic  representation 


f  -2 

X—  0 

a  m„ 


0 


0  Yj—  0 


«  ma 


0 


0 


+ 


fV  V  0A 

xx  xz 

V„  V„  0 


xz  zz 


OOF 


yy  J 


«  ™aj 


f  F  \ 

E 

0 

0' 

f  F  ^ 

= 

0 

E 

0 

KFyy) 

0 

0 

E 

J^yy) 

where  the  NACTs  have  been  eliminated  and  the  diabatic  wavefunctions  are  coupled  by 
the  diabatic  PES. 


The  Asymptotic  Representation 

To  facilitate  wavepacket  propagation  as  required  by  the  channel-packet  method 
(described  in  a  later  section),  the  Hamiltonian  can  be  represented  in  an  asymptotic  basis. 
Before  introducing  this  basis  it  is  necessary  to  establish  a  new  coordinate  system. 

Body  Fixed  and  Space  Fixed  Coordinates 

The  coordinates  used  by  Weeks  et  al  [16]  are  shown  in  Figure  2. 


(27) 
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Figure  2.  Space  and  Body  Fixed  Coordinate  System  for  B+FI2 


The  B  +  H2  system  is  described  using  the  Jacobi  coordinates:  R  is  the  distance 
from  the  boron  atom  to  the  center  of  mass  of  the  hydrogen  molecule;  r  is  the  hydrogen 
bond  length;  9  is  the  angle  fonned  between  the  hydrogen  axis  and  the  line  connecting  the 
boron  atom  to  the  hydrogen  center  of  mass.  The  space  fixed  axes  are  labeled  X,  Y,  Z 
while  the  body  fixed  axes  are  labeled  x,  y,z .  Both  have  their  origin  at  the  B  +  H2  center 

of  mass,  which  in  the  figure  is  displaced  away  from  the  boron  atom  for  visualization 
purposes.  The  body  fixed  z  axis  contains  both  the  boron  atom  and  the  hydrogen  molecule 
center  of  mass.  The  azimuthal  angle  <f>  is  fonned  between  the  body  fixed  x  axis  and  the 
projection  of  the  hydrogen  molecule’s  axis  onto  the  xy  plane,  and  is  omitted  from  the 
figure  for  clarity.  The  body  fixed  y  axis  is  constrained  to  lie  in  the  XY  plane.  The  body 
fixed  axis  system  is  oriented  in  the  space  fixed  coordinate  axes  using  the  Euler  angles  a 
and  /?.  The  angle  a  is  defined  between  the  projection  of  the  z  axis  onto  the  XY  plane  and 
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the  X  axis.  fl  is  the  polar  angle  fonned  between  the  body  fixed  z  axis  and  the  space  fixed 
Z  axis. 


Using  these  space  fixed  and  body  fixed  coordinates,  the  Hamiltonian  may  be 
rewritten  in  the  center  of  mass  frame  [17]  as 


H 


_  Pr 


J 


Is 


^  Mb,  h,  2//#  2pHr  2//  B,H 


2+H?+H*  +  Vel+Hso+Eoff 


The  first  four  terms  in  eq.  (28)  represent  the  nuclear  kinetic  energy  operators  using  the 


reduced  mass  of  the  nuclei  pB  H  = 


2  mBmH 
mB  +  2  mH 


,  the  reduced  mass  of  the  hydrogen 


molecule  juHi  =  — — ,  and  the  momentum  operators  for  their  respective  coordinates,  pR 


and  pr .  The  angular  momentum  operator  j  describes  the  II 2  rigid  rotor,  while  L 
corresponds  to  the  rigid  rotor  composed  of  the  boron  atom  and  the  hydrogen  molecule 
center  of  mass.  H^2  and  represent  the  kinetic  energy  and  coulomb  potential  of  the 

electrons  in  the  molecule  and  the  atom,  while  Vel  represents  the  coulomb  potential 
between  the  boron  atom  and  the  hydrogen  molecule.  Hso  is  the  spin-orbit  Hamiltonian 
for  the  boron  atom,  and  Eoff  is  an  arbitrary  energy  offset. 


A  complete  basis  for  this  system  is  given  by 


J  ,l,sja  ,r,R)  where  the 
MP  k  CD  / 


coupling  scheme  is  identified  as  “case  1A”  by  Dubemet  and  Hutson  [17]. 


(28) 
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Figure  3.  Angular  Momentum  Coupling  Scheme 


The  electronic  orbital  angular  momentum  /  and  spin  s  of  the  boron’s  unpaired 
outer  electron  couple  to  give  ja  with  projection  <x>  onto  the  body  fixed  z  axis.  The 
angular  momentum  of  the  H2  rigid  rotor  is  j  with  a  projection  k  onto  the  body  fixed  z 
axis,  and  there  is  assumed  to  be  no  coupling  between  j  and  ja .  The  total  angular 
momentum  is  J  =  L  +  j  +  ja  ,  which  in  the  centrifugal  sudden  approximation  [18]  has  a 

projection  onto  the  body  fixed  z  axis  of  P  -k  +  co  and  a  projection  onto  the  space  fixed  z 
axis  of  M. 


This  basis 


J  j  j  a  \ 

,  ,  /,  s,  a  ,r,R)  has  an  infinite  number  of  states,  but  may  be 
MP  k  co 


truncated  to  include  only  the  range  of  states  accessible  as  determined  by  the  kinetic 
energy  of  the  B  +  H2  system.  The  electronic  basis  has  been  truncated  to  only  include  the 


boron  2  P  states,  which  gives  1  =  1,  s-  1/2  which  limits  ja  to  the  values  1  /  2, 3  /  2  . 
Additionally  the  hydrogen  molecule  is  restricted  to  the  ground  electronic  state.  This 


15 


work  considers  only  a  total  angular  momentum  J  —  1  /  2  .  A  more  concise  notation  for 
the  basis  is  obtained  by  including  only  those  labels  which  are  not  constant. 


J  J 
M  P’k 


1 A 


co 


Representing  the  Hamiltonian 

The  Hamiltonian  in  eq.  (28)  is  now  represented  in  this  basis.  The  details  behind 
the  results  presented  here  are  covered  in  more  detail  by  Weeks  [16]  and  again  by  Garvin 
[1].  The  matrix  elements  of  the  first  three  terms  in  eq.  (28)  are  diagonal  in  this  basis. 
They  are 


where  5X,X 
diagonal. 


,r',R' 

*2 

Pr 

\k'  co' 

2Pb,h2 

j  Ja  D\  -h2  1  DS 
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h  i 

J  ,a  ,r',R 
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~  2 
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JJ-,r,R)  =  dL  liLrf, 

k  co  j  2juh  r  dr~ 


I  j'  j  ’ 

j  Ja  l  nl 
,  ?  ?  1  ?  -*v 


\k'  co' 


2Ph/~ 


7  Ja  r,p\=  h  J  S 


k  co 


2PhJ 


2  xx 


Sp'pS'hi S'Jm  8{r-  r ')  5{R  -  R ')  .  The  remainder  of  the  terms  are  not 


(29) 


(30) 


(31) 


The  fourth  term  is  obtained  by  solving  for  L  in  terms  of  the  total  angular 
momentum  J ,  then  expanding  using  raising  and  lowering  operators. 


■Jl  +j2  -J+h_  -J  Ja,  -  2-]Ja  +jja  +  2jJa . 


(32) 


Each  tenn  in  eq.  (32)  is  evaluated  to  give  the  full  expression 
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J  Ja  r » 

9  9  '  9  -*v 

\A: 1  ft)1 
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J  Ja 
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,r9R)  = 


I  j'  j  1 
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-  , 9  9  '  9  ±V 

\k  co 


(■ J-J-Ja ) 


2  Mb.hR2 


J  Ja 


9  9  7  9 


k  co 


r.R 


2Mb.hR 2 


{j(j  +  i)  +  j(j  +  i)  +  ja  (ja+\)-2(k  +  a>)  +2 ka>}sx,x 


{{j(j  +  ^)~P(P  ml))(y  (y  +  \)-k(k  ml))}' 2  SJMJ,MS r  Sk,kxA8 jd  j8 m,a8{R  -R') 

{( j(j  +  l)-P(Pml))(ya(ya+l)-®(®ml))}1/2^' fuSrAA'j.s»™s(R~R') 

{( ^  (^  + 1)  -  A:  (A:  ±  1))  (ya  + 1)  -  « (fi»  ml))}1/:!  v.  ^  0 


The  second  and  third  terms  in  eq.  (33)  are  neglected  under  the  centrifugal  sudden 
approximation. 


The  spin-orbit  Hamiltonian  is  given  by  Hso  =  t^l  -s ,  where  C,  is  treated  as  a 

constant  under  the  pure  precession  approximation  [6],  and  when  represented  in  the 
angular  momentum  basis  is  given  by 


H 

\k'  (O' 

SO 

J  J , 


k-a,-r'R)=2  Ub.+O-^  +  O-sb  +  OK.  (34) 


The  remaining  terms  are  Hd2  +  Vel  +  HR  +  Eoff .  These  are  represented  in  the 
asymptotic  electronic  basis  |Z,A)  =  |  'Z^®|/,  A)  consisting  of  the  ground  state  of  the 
hydrogen  molecule  electronic  Hamiltonian  |  *Z*^  and  the  three  2  P  states  of  the  boron 
electronic  Hamiltonian,  given  by  |/,A)  with  possible  values  /  =  1,  A  =  -1,0,1 .  In  the 

asymptotic  limit  the  term  (Z ',  A  |  Hd  |  Z,  A^  is  constant,  and  Eoff  =  -  ^Z ',  A  |  Hd  |  Z,  A^  is 
chosen  to  eliminate  these  terms. 


(33) 
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/  j  '  J  a  ' 

The  representation  of  the  remaining  terms  {  ,  a  ,r\R' 

\k'  co' 

accomplished  using  the  numerically  detennined  diabatic  electronic  potential  energy 
surfaces  from  equation  (27).  These  are  cast  in  a  new  form  by  expanding  in  terms  of 
reduced  rotation  matrix  elements  which  are  related  to  the  associated  Legendre 

polynomials  Plm  (cos  6)  by  dlm0  ( 6 )  =  [(/  - m)\l  (/  +  /«)!]'  "  P‘n  (cos 9 ) .  This  expansion 

takes  the  fonn 


H"2  +K, 


J  Ja 
k’  co 


j\R)  is 


r”M,0)=£  v^(r,R)d%(e) 
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KPR-e)=£  P(r,R)d^(ff)-£  V^(r,R)di(6) 


Ar=  0 


A**  — 2 


The  expansion  coefficients  V'j ,  Vx,r ,  V/  ,  fy  in  equation  (35)  are  obtained  by 

perfonning  a  multipole  expansion  of  (  S', A’  |  Vel  +  // J2  |  X,  A  ) .  Integrating  over 

molecular  and  atomic  electronic  coordinates,  with  the  exception  of  the  polar  and 
azimuthal  angles  of  boron’s  outer  2 p  electron,  gives  the  matrix  elements  [6] 

CL-mH^.a'IA+A?  |z,a> 

=  (r,A'|£,-,  py  |  S,A)  +  (I',A'|i/J-|  S,A) 
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(35) 


(36) 
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Here  the  terms  indicated  by 


are  the  Wigner  3j  symbols  and  (0,<f>)  are 


V-  •  -J 

modified  spherical  hannonics.  Collisional  energies  considered  here  will  not  allow  an  H2 
electronic  transition  so  may  be  removed.  The  value  Aa  is  restricted  to  Aa  =0,1,2  by 
the  triangle  inequality  |l  -  Aa  |  <  1  <  |l  +  | .  This  reduces  equation  (36)  to 
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(37) 


A  further  restriction  //  =  A  -  A '  is  obtained  by  considering  the  bottom  row  of  the  Wigner 
3j  symbols,  where  it  is  required  that  -A fi  —  A  =  0 .  Defining  a  new  form  of  expansion 

coefficients  VA  -  ^  C';  (0,O)VA  A  M(r,R )  the  ju  matrix  elements  are 

x. 
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Transfonning  to  a  Cartesian  basis 


X>  =  ^(-I1>  +  |-I))’k) 


(l1>  +  |-1».k>  =  l°> 


gives 
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(39) 


VD(r,R,0)  = 
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Each  matrix  element  in  equation  (39)  is  compared  to  the  expansion  of  the 
numerical  diabatic  potential  energy  surfaces  in  equation  (35),  for  which  the  fit  process 
has  given  expansion  coefficients  Vzir ,  ,  Vfr ,  Vy  as  a  function  of  r  and  R.  These  are 


related  to  the  expansion  coefficients  VX  /  fi  (r,  Rj  by 


^(r.R)=\[V:f(r.R)  +  2V^(r,R)]-VHi(r)Sv 
VX2M’R)  =  l[Kr  (r,R) 


(40) 


H , 


Having  obtained  the  expansion  coefficients  the  interaction  electronic  potential  Ve/  +  Hel 
is  represented  in  the  angular  momentum  basis.  This  is  given  by  [16] 

I  j  J'a 
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H7  +K, 
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[ja'kja\ 

'  h  K  ja' 

y-k'  /./  ky 

v0  0  0 J 

1  /  J  /  J 

v  co  —fi  —co'y 

,0  0  0, 

S"„Sp.r5irVS(R'-R) 


(41) 


where  terms  of  the  form 


(L  )  [L  1 

are  3  -  /  symbols  and  those  of  the  fonn  <  }  are  6-  / 

h  J  M 


symbols. 
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The  Effective  Diabatic  Potential  Energy  Surfaces  for  B+H2 

Taken  together,  eqs.  (31),  (33),  (34)  and  (41)  give  the  effective  diabatic  potential 
energy  matrix  elements 


VeDff(i\R)  .  =(j  ,Ja  ,r\R 

eff\  )yy  \k'  Q)' 


-  +  - 


z2 


2/uHr-  2  fi  R 


2+H:2+vel+Hsi 


J  Ja 


?  n '  n 


k  a> 


r,R 


(42) 


where  the  label  y  represents  appropriate  values  of  J,  M,  P,j,  k,  ja  and  co.  These 

represent  surfaces  for  which  the  time-dependent  Schrodinger  equation  can  be  solved. 
The  surfaces  are  two  dimensional  (for  r  and  R)  and  each  is  labeled  by  the  choice  of  y . 
The  diabatic  effective  potential  energy  surfaces  used  in  this  work  were  calculated  by 
Garvin  [1], 


The  full  Hamiltonian  is  obtained  by  adding  the  nuclear  kinetic  energy  tenns  from 
eqs.  (29)  and  (30)  to  the  effective  potential  in  eq.  (42).  Reduced  wavefunctions 
0(r,R,t)  =  rRW  ( r,R,t )  are  introduced  to  simplify  the  radial  form  of  the  kinetic  energy 
operators,  giving  the  time-dependent  Schrodinger  equation  as 
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0  J 
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(43) 
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Eq.  (43)  shows  that  each  wavepacket  is  dependent  on  more  than  one  effective  diabatic 
PES  Vej  as  it  evolves  in  time. 

Structure  of  the  Effective  Diabatic  Potential  Energy  Matrix 

The  effective  diabatic  potential  energy  matrix  V^f  (r,  R)  ,  is  block  diagonal  as 

governed  by  the  Kronecker  deltas  in  eqs.  (31),  (33),  (34)  and  (41).  These  give  the  matrix 
a  structure  of  infinite  dimensional  blocks  arranged  in  a  hierarchy.  The  largest  are  labeled 
by  their  value  of  total  rotational  energy  J.  Each  J  block  is  subdivided  into  identical 
blocks  labled  by  M.  This  study  is  restricted  to  J  —  1  /  2  .  This  block  is  shown  in  Figure  4. 


V 

Figure  4.  Structure  of  the  J  =  1/2  Block  of  the  Effective  Potential  Energy  Matrix 


Since  the  M  blocks  are  identical,  only  M  =  1/2  is  considered.  Under  the  centrifugal 
sudden  approximation  [  1 8]  AP  =  0  ,  so  we  may  consider  only  the  P  =  1/2  block. 
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The  J  =  \!  2,M  =  H  2,P  -  2  block  is  further  divided  into  orthohydrogen  and 


parahydrogen  blocks.  In  this  work  only  transitions  from  the  ground  rotational  state  j  —  0 
are  considered,  so  only  the  even  valued  j  states  corresponding  to  the  parahydrogen  block 
are  considered.  The  requirement  k  +  (o  =  P  =  \/2  along  with  ja  =  T,l  restricts  the 
number  of  states  available  for  each  value  if  j.  Specifically  they  are 


7  =  0 
7  =  2,4,... 


j  Ja\_ 

°  *\ 

k  co  / 

0  V 

0 
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7  Ja 
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1  ~h 
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1  ~h 


j  1\ 
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0  h 
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-l 


The  total  number  of  these  states  is  infinite  but  may  be  truncated  to  only  include  those 
states  which  are  energetically  accessible  as  determined  by  the  total  energy  of  the 
wavepacket  being  considered.  Truncating  to  only  include  values  up  to  ymax  results  in  the 
effective  potential  energy  matix  having  a  dimension  given  by 


(44) 


n  =  2  +  3ymax  (45) 

This  n  is  also  referred  to  as  the  basis  size,  since  the  effect  of  truncating  is  to  restrict  the 
number  of  basis  states  in  which  we  represent  the  Hamiltonian  to  n. 

As  seen  in  Eq.  (43)  the  off-diagonal  terms  of  V°ff  ( r,R )  couple  the  evolving 

wavepackets,  allowing  for  transitions  between  the  various  surfaces.  These  transitions 
represent  a  wavepacket  in  an  initial  channel  y  changing  to  a  final  channel  ;/ '  as  a  result 
of  scattering.  The  tenn  channel  refers  to  a  specific  arrangement  of  B  +  H2  together  with  a 
choice  of  quantum  numbers.  The  process  for  calculating  the  probability  of  making  these 
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transitions  is  the  subject  of  subsequent  sections  on  scattering  and  the  channel  packet 
method. 


Nonadiabatic  behavior  is  most  frequently  associated  with  coupling  between 
electronic  and  nuclear  degrees  of  freedom.  Representing  the  Hamiltonian  in  the  angular 


momentum  basis 


J  J 
M  P’k 


CO 


has  resulted  in  an  effective  potential  energy 


matrix  VeZ  which  contains  off-diagonal  tenns.  This  nonadiabaticity  is  the  result  of 

coupling  between  rotational  and  vibrational  degrees  of  freedom  as  well  as  coupling 
between  electronic  and  nuclear  degrees  of  freedom. 


The  Split-Operator  Propagator 

The  time-dependent  Schrodinger  eq.  (2)  gives  the  solution  to  the  initial-value 
problem  [19] 

itH 

i//(x,t)  =  e  h  ^(*,0)  =  Ui//(x,0)  (46) 

In  many  cases  this  general  form  of  the  evolution  operator  U  must  be  approximated.  The 

split-operator  method  is  one  technique  for  calculating  U ,  thereby  giving  the  time 
evolution  of  an  initial  state  based  on  knowledge  of  the  system’s  Hamiltonian. 


The  Split-Operator  Approximation 

The  split-operator  method  uses  the  approximation  [20] 
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(47) 
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where  0(  A/  )  represents  terms  that  are  third-order  or  higher  in  At . 


The  equivalence  of  this  approximation  up  to  third-order  or  higher  terms  may  be 
shown  by  considering  Taylor  series  expansions: 
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As  shown  in  eq.  (47),  the  split-operator  approximation  ignores  all  tenns  that  are  third- 
order  or  higher  in  At .  Therefore  the  term  0(  At3 )  is  the  error  associated  with  this 

approximation,  which  can  be  made  arbitrarily  low  by  selecting  a  sufficiently  small  time 
step  At. 
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Diagonal  Representations 

The  motivation  for  this  approximation  is  the  relative  simplicity  gained  by 
representing  each  operator  in  a  basis  in  which  it  is  diagonal  [20], 

When  the  effective  potential  V  is  not  diagonal,  as  is  seen  for  the  effective  diabatic 
PESs  of  B  +  Ht,  the  matrix  representation  of  V  can  be  diagonalized  using  a  unitary 

transformation  [21]  into  the  adiabatic  representation.  This  transformation  will  depend  on 
the  coordinates  being  propagated  over  -  in  the  case  of  the  diabatic  effective  PESs  of 
B  +  H2,  the  transfonnation  may  be  represented  by  UAD  (r,  R) ,  and  the  potential  is 
diagonalized  at  each  point  on  the  grid  by 

(48) 

The  wavepacket  is  also  transformed  from  the  diabatic  to  adiabatic  representation  by 

(49) 

If  the  matrix  representation  of  V  and  f  are  diagonal  in  the  appropriate 
representation,  so  will  the  matrix  representation  of  the  exponential  of  those  operators, 

— —VAt  --TAt 

e  2h  and  e  h  .  Recall  that  an  exponential  with  a  matrix  in  the  argument  is  defined  by 
its  Taylor  series  expansion.  Successive  powers  of  a  diagonal  matrix  will  be  diagonal,  and 
the  series  may  be  reduced  to  a  diagonal  matrix  whose  elements  are  the  exponent  of  the 

original  matrix.  For  example,  if  V is  the  matrix  representation  of  V  in  a  basis  where  it  is 
diagonal, 

%  0  0  N 
V=  0  F22  0 

l  0  0  Ks, 
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(50) 


Changing  Representations:  The  Fourier  Transform 

Assuming  the  potential  energy  operator  V  and  the  kinetic  energy  operator  T  are 
functions  only  of  the  spacial  coordinates  and  momentum  respectively,  V  is  diagonal  in 
the  coordinate  adiabatic  representation,  while  T  is  diagonal  in  the  momentum  diabatic 
representation. 

To  show  this,  consider  the  matrix  element  of  a  potential  which  depends  only  on 
the  coordinate  x,  (x'|F|x^ .  Since  V  is  only  a  function  ofx,  we  may  expand  it  as  a  series 
in  terms  of  x  . 

(jc'|F|jc)  =  (jc'|C0  +Clx  +  C2x2  + ...  |  x'j 

=  (x’|C0  Ix^  +  ^x'ICjxlx^  +  ^x'lQx2 1 x'j  + ... 

Since  the  x)  are  eigenstates  of  the  x  operator,  this  is 

(x'|K|x)  =  C0  (x'|x)  +  C1x/  (x'|x)  +  C2x2  (x'|x)  +  ... 

(x’|F|x)  =  (C0+C1x  +  C2x2  +...)S(x-x')  (51) 


Since  the  eigenstates  of  the  position  operator  are  orthogonal,  the  matrix  element  is 
nonzero  only  if  it  is  on  the  diagonal.  The  same  may  be  shown  for  (k '  |  T  |  k'j . 
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Similar  to  how  the  unitary  transfonnation  UAD  is  used  to  change  between  the 

diabatic  and  adiabatic  representations,  there  is  a  transformation  to  move  between  the 
coordinate  and  momentum  representations  -  the  Fourier  transfonn.  To  understand  how 
this  is  done,  consider  the  coordinate  representation  of  our  time-evolution  operator  on  an 
initial  state: 

/  I  -iV AiV2h  -iTA//h  -iVAt/2h\,,\ 

\x\e  e  e  \v) 

Using  the  completeness  relation  j* 1  x)  (x  |  dx  =  J  |  k}  (k  |  dk  =  I ,  where  I  is  the  identity 
operator,  we  may  rewrite  this  as 


-;TA//2h  -iTAtlh  -iVAtllh 

A  C-  t-  C/ 
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=  ( 
=  { 


-iVAtHh  ff-iTAt/h 


x\e  "  ~‘IIe 

iVAt/lh 


lie1 


FA/72  h 


7>) 


x  e 


J|x  ^(x  k^(k  e  7A,/hJ|A:^A:|J|x  ^(x  e  '*  A'/2h  J|x  ^x  y/)dx  dk  dk  dx  dx 


This  may  be  rearranged  by  moving  terms  that  are  not  functions  of  the  respective 
integration  variables. 

(x|e-'fAj/2V'fA,/V^/2h|  </) 

=  j*(x|e~'’  A'/2h  |x  |A:)|(A:  ,h \k^j  J(^|x  ^x  |u'(  A'/2h  |x  ^x  | y/^ dx  dk dkdx  dx 

Eqs.  (5 1)  and  (50)  have  shown  that  the  terms  (x  |  e“,KA'/2h  |x)  and  (k  |  e“,rA'/h  |  k)  are 
diagonal,  that  is,  they  contain  the  delta  functions  £ (x - x  )  respectively.  This 

eliminates  three  of  the  integrals  giving 

(x\e-l™l2he-i™l\-ifMI2'n\i//) 

=  (x  |  e~l * A'/2h  |  x) |  ^x  |  k  ^(k  |  e~lfA>/  h  |  k'j  J (k  |  x  ^x  |  e~'1  A,/2h  |  x  'j  ^x  |  y/^dkdx 
Recalling  (xt  |  x  .  ^  =  (kj  |  U  )  =  djJ,  this  is 
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(x 

(x 


1  -iV&t/2h  -iTAtlh  -iFAf/2h  1  \  _ 

1  W  W  w  1  VIS  / 

I e-^M,2h  | x)  J (x \k' }  dx  {k  | e-ifA,,h  | k)  j 

[^kjx^Jk^x  e  ,fA'/2h  x^x  | 

(52) 


Here  we  can  pause  and  identify  these  terms,  (x  |  e  lVAt,2h  |  x)  is  the  coordinate 
representation  of  the  first  part  of  the  time-evolution  operator,  (k  |  e  ,rA'/h  |  is  the 
momentum  representation  of  the  second  part.  The  remaining  terms,  J  ^x  |  k ^  dk  and 
^{k  |x ^  dx  ,  represent  the  transformation  between  the  two  representations.  In  fact, 
recognizing  (x\k  'j  as  the  coordinate  representation  of  the  momentum  eigenstate,  which 

are  plane  waves,  we  can  rewrite  this  as  — f  elkxdk  ,  meaning  we  can  cast  this  as  an 

V2  J 

Inverse  Fourier  Transform,  and  J^k|x  ^  dx  ,  its  Hermitian  conjugate,  as  the  Fourier 
Transform. 

This  is  an  important  result  in  terms  of  ease  of  computation.  In  practice, 
transformation  between  the  two  representations  is  typically  accomplished  using  a  Fast 
Fourier  Transform  (FFT)  algorithm.  On  a  grid  of  size  N,  the  FFT  requires  0(  N  log2  N  j 
operations  [22],  making  this  propagation  scheme  viable  for  relatively  large  grids. 

The  Split-Operator  Calculation 

Putting  the  transformations  described  in  eqs.  (48)  and  (52)  together,  if  given  a 
wavepacket  in  the  diabatic  coordinate  representation  at  time  t  -  0,  (/>r>  (r,R,  0) ,  first 
transform  this  to  the  adiabatic  representation 

Uj{r,R)tD{r,R,  0) 
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Next  apply  the  first  of  the  potential  operators  e  ' 

e-ivAAtnhUj(r,R)r(r,R,  0) 

where  VA  is  the  diagonal  adiabatic  coordinate  representation  of  the  potential.  This  is 
transformed  back  into  the  diabatic  representation  and  then  moved  into  the  momentum 
representation  via  the  Fourier  transform,  given  by  the  operator  F 

FUADe-'‘^UAD'(r,R)r(r,R,0) 

This  is  multiplied  by  the  exponential  kinetic  energy  operator  e“,r°A'/h  where  TD  is 
diagonal  in  the  diabatic  momentum  representation,  then  moved  back  to  the  coordinate 
representation  via  the  inverse  Fourier  transfonn  F 

F-'e-iTD^FUADe-iMUj  {r,R)fD  (r,R,0) 


The  wavepacket  is  again  transfonned  to  the  adiabatic  representation  and  multiplied  by  the 
adiabatic  representation  of  the  second  of  the  exponential  potential  operators 

JF-xe-iT^FU ADeivA^lj  J  (r,R)f  (; r,R,0 ) 


Finally  this  is  transformed  to  the  diabatic  representation  to  yield  the  wavepacket  at  t  -  At 


r  (r,R,  At)  =  UADe-ivA  J  F^e^FU ADe~iMU J  (r,R)f  (r,R,  0)  (53) 


The  choice  to  begin  and  end  with  the  wavepacket  in  the  diabatic  representation  is 
arbitrary.  Propagating  in  the  adiabatic  representation  eliminates  the  first  and  last 
transformations  from/to  the  diabatic  representation. 


4>A  ( r,R,At )  =  e-ivAM,2hUAJF-le-iTAt,hFUADe-iV*A,l2h  ( r,R)</>A  (r,R,  0) 


iVAMI2U  , 


(54) 
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Scattering  Theory 


The  split-operator  propagator  allows  us  to  take  a  wavepacket  in  an  initial  known 
state  and  propagate  it  in  time.  This  tool  may  be  applied  to  examine  the  behavior  of  a 
wavepacket  used  to  describe  molecular  dynamics  as  governed  by  a  suitable  potential. 
Tracking  possible  changes  of  state  in  the  wavepacket  as  a  result  of  these  dynamics  is  the 
domain  of  scattering  theory. 

Asymptotic  Reactants  and  Products 

Experimentally,  the  actual  quantum  scattering  event  takes  place  on  a  scale  that  is 
unobservable.  As  a  result,  instead  of  discussing  the  interaction  event  scattering  theory 
relates  widely  separated  reactants  located  in  the  asymptotic  limit  before  a  collision  to 
products  located  in  the  asymptotic  limit  after  the  collision.  A  reactant  wavepacket  in  the 
asymptotic  limit  before  coming  into  contact  with  the  potential  will  be  called  \y/in)  ,  while 

an  asymptotic  product  wavepacket  is  called  \y/out) . 

The  Moller  Operators 

There  are  several  definitions  which  aid  in  this  process  of  relating  initial  states  to 
asymptotic  wavepackets.  The  full  Hamiltonian  of  the  system  can  be  broken  into  two 
parts,  H  =  H0  +  V  ,  where  V  contains  all  terms  describing  the  interaction  between 
reactants  or  products  [2],  The  tenn  H0  is  the  asymptotic  Hamiltonian  describing  the 

system  in  the  asymptotic  limit  when  V  — »  0 ,  that  is,  where  there  is  no  interaction  between 
reactants  and  products.  The  full  and  asymptotic  Hamiltonians  are  used  to  define  a  set  of 
isometric  operators,  known  as  Moller  operators  [23]: 

Q+  =  lim  e'^V'A/.'h 

t->-  OO 


(55) 
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(56) 


Cl 


These  operators  are  made  up  of  two  time-evolution  operators.  The  result  of  Q+  acting  on 
a  state  \i//)  is  to  first  propagate  it  backwards  in  time  to  t  =  -go  under  the  asymptotic 
Hamiltonian  H0 ,  then  to  propagate  it  forwards  in  time  to  t  =  0  using  the  full  Hamiltonian 
H  .  Similarly  the  effect  of  Cl  on  the  state  |  (ft)  is  to  propagate  it  forwards  in  time  to 
t  =  co  under  H0 ,  then  to  propagate  it  backwards  in  time  to  t  -  0  using  H  . 


The  Intertwining  Relation 

These  two  Hamiltonians  and  the  Moller  operators  satisfy  the  intertwining  relation 

[23]: 

hCi+  =  Cl+H0 


To  see  this,  consider 


iHz/h/-\  iHr/h 

e  Ll.=e 


1  • _  mtm 

lim  e  e 

z->mx) 


iHt!  h  -iH0t/h 


=  lim 

/^m» 


giH  (r+t)/h  g-iH0t/h 


_  eiH(T+t)/he-iH0(T+t)lhgiH„Tlh 


=  fl+e 


iHnT/h 


Differentiating  with  respect  to  r  and  then  setting  r  =  0  yields  eq.  (57). 

The  Moller  States 

The  new  states  obtained  by  applying  the  Moller  operators  to  a  state  at  t  -  0  are 
called  the  Moller  states: 


(57) 
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(58) 


k+)=6+  \w) 

k-)=^-W 

The  M0ller  operators  give  the  actual  states  at  t  -  0  that  would  evolve  from  (or  to)  the 
asymptote  represented  by  the  state  they  operate  on  [23].  Each  state  \y/)  at  t  -  0  is  related 

to  an  in  or  out  asymptote  \y/in/out)  by  the  Moller  operators: 

k)  =  «>*}  =  «>„,)  (59) 

The  Scattering  Matrix 

The  Moller  operators  are  isometric,  so  eq.  (59)  may  be  inverted. 

WoUt)=^-W)=^AWin)  (60) 

The  scattering  operator  is  defined  as 

S  =  iYo+  (61) 

The  scattering  operator  relates  asymptotic  reactant  states  to  asymptotic  product  states. 

\v0u,)  =  S\Vin)  (62) 

If  the  Moller  states  are  obtained  from  arbitrary  product  and  reactant  states  at  t  -  0  as 
given  by  eq.  (58),  the  overlap  of  the  Moller  states  represents  the  probability  amplitude 
that  the  reactant  state  \y/}  will  scatter  into  the  product  state  \(f)  [23].  This  overlap  is 
given  by 

k_k+)  (63) 


Recalling  the  definition  of  the  Moller  states,  this  may  be  written  as 

{y/_\y/+)  =  {(p\Cl'_Cl+\y/) 

=  {(j)\S\y/) 


(64) 
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So  the  matrix  element  (<j)\S\y/^ gives  the  probability  amplitude  that  \i//)  will  scatter  to 


IA 


Representing  the  States 

A  common  choice  for  representing  the  states  | y/^  and  \(/>)  is  the  basis  of 


eigenstates  of  the  asymptotic  Hamiltonian  H0 ,  which  we  will  label  as  ky ,  y'j  where  y 

identifies  the  quantum  numbers  associated  with  all  internal  degrees  of  freedom.  These 
eigenstates  have  eigenvalues  given  by 


(65) 


h  X 

where - is  the  relative  kinetic  energy  and  E  is  the  energy  associated  with  internal 

2m 

degrees  of  freedom,  and  E  is  the  total  energy. 

Another  set  of  states  that  will  later  prove  useful  is  obtained  by  operating  on  this 
original  basis  with  the  Moller  operators: 


(66) 


These  new  states  are  eigenstates  of  the  full  Hamiltonian  H  .  From  the  definition  in  (66) 

\kr>r) 
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r  +E„ 
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(67) 
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Also  using  the  intertwining  relation 


Ci±H0  kr,y)  -  HCi±  kr,y 


±  r’ 


n±H0  kr,y)  =  H  kY,y± 


(68) 


Taken  together,  (67)  and  (68)  show 

H\kr±,y)- 


f  h\2  ' 

~W^  +  Er 
v  2m  y 


kr,y±) 


(69) 


which  shows  k±,  y)  are  eigenstates  of  the  full  Hamiltonian. 


The  product  and  reactant  states  are  represented  in  the  k  ,  y)  basis  as 


\v)  =  J  *1+(kr)\kr,y}dkr 


(70) 


| <f)=\  n~(kr)\kr,y)dk: 


Using  this  we  can  write  the  Moller  states  as 

k±)  =  n±k> 


00 

=  Q±  J  n±(ky)\k7,y)dky 

-oo 

00 

=  J  T(kr)n±\kr,y)dkr 

-00 


W±)=  J  n±(ky)\kr,y±)dkr 


(71) 


(72) 
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Note  the  77's  are  the  same  in  equation  (72)  as  in  (70)  and  (71)  -  that  is,  the  expansion 
coefficients  for  the  Moller  states  in  the  basis  of  eigenstates  of  the  full  Hamiltonian  are  the 
same  as  those  of  the  asymptotic  product  and  reactant  states  in  the  basis  of  eigenstates  of 
the  asymptotic  Hamiltonian. 

These  basis  states  satisfy  the  orthogonality  relation  [2] 


(kr.,r'-\kr,r+) 


S(E'-E)S^k> 


(73) 


where  S[:[  is  the  on-shell  S-matrix  element.  This  element  gives  the  probability 


amplitude  to  scatter  from  an  initial  state  with  momentum  ky  and  internal  state  y  to  a  state 


„  h  V  ^ 

with  momentum  k  ,  and  internal  state  y ' .  Rewrite  this  using  E  - - h  E  giving 

2m 


This  may  be  recast  using  two  properties  of  the  Dirac  delta  function: 

<5(ax)  =  ald  (x) 


(75) 


S( 
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(76) 


Rewriting  (74)  with  (75)  and  (76)  gives 
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1/2 

(kfy-\kr,r+)=-^sift\s(kr-k,)+s(k,+kr)) 


(77) 


These  relations  are  instrumental  in  developing  a  technique  for  calculating  scattering 
matrix  elements  based  on  the  correlation  function  between  evolving  reactant  and  product 
Moller  states.  This  is  known  as  the  channel  packet  method. 


The  Channel  Packet  Method 

The  final  piece,  calculating  the  S-matrix  elements,  has  been  foreshadowed  by  the 
orthogonality  relationship  expressed  in  (77).  This  technique,  known  as  the  Channel 
Packet  Method  (CPM),  is  laid  out  by  Weeks  and  Tannor  [2], 

Begin  by  considering  the  Fourier  transform  of  the  Moller  state  \ys+),  written  in 
time-dependent  fonn: 


|4(£))=  je'E'"'\wJt)}dt=  J 


iEtlh  -iHl/h 

e  e 


y/+)dt 


(78) 


Using  the  expansion  given  by  (72), 
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The  second  part  of  (79)  represents  a  delta  function 


(80) 


Using  (80),  along  with  (75)  and  (76)  rewrite  (79)  as 


K(*))=  j  '?*(M|*,.r+)(tf(%-id  +  'J(tr'+tr)))  dkr 


\A+(E))=mf 


'(-ky)\-ky,r+)) 


(81) 


Next  evaluate  the  scalar  product  (i//  At  (Efj  ■  Using  (72)  to  expand  the  Moller  state  this 


is 


(W- 1 A  (E)}  =  j  if*  (kr )  (kr;  7 1  ^  [rf  (+kr)  |  +ky ,/+)  +  ri+(-kr )  |  -kr ,  y  +))  dkf 


=  ^\\W(kr)^+^+k^{kr’y'~\+ky’y+)  +  Tl*(kr)ri+  k ~kr )  (kr ’  ^ u  I  > r+))dky  (82) 


Using  the  orthogonality  relation  from  (77)  this  is 
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(83) 


{v- 1 A  (E))  =  [>r‘  (-*,)'?* (+kr)S%Mj 

+’l"(+kr)’?'(+WC,.,, 

+A'{-K)1  \-K)S%,-lr 
+n-(+kr)„'(-kr)s:(^r\ 

Recall  that  the  expansion  coefficients  77*  {ky )  simply  describe  the  momentum 

representation  of  the  initial  wavepacket  we  wish  to  propagate.  Hence  we  are  free  to 
impose  certain  restrictions  which  will  simplify  eq  (83).  For  example,  we  can  choose  to 
select  an  asymptotic  product  state  which  only  has  positive  momentum  -  that  is,  we 
impose  the  restriction 


(84) 


It  is  worth  noting  that  apart  from  being  mathematically  expedient,  choosing  an  incoming 
wavepacket  with  entirely  positive  momentum  is  practical  [2],  as  instead  of  spreading  out 
in  all  directions  all  of  the  wavepacket  will  propagate  towards  the  interaction  region  we 
wish  to  examine.  Similarly,  reactant  wavepackets  are  chosen  to  have  only  negative 
momentum: 

?l+(+kr)  =  0  (85) 

Imposing  the  restriction  from  (84)  onto  (83)  results  in  all  but  the  second  term  on 
the  right  side  being  zero: 

<f  K)  =  4  (E))  =  Kb*  (+kr)S%„  (86) 

hvrHn'l] 

Alternately,  by  choosing  appropriate  combinations  of  the  asymptotic  product  and/or 
reactant  states  containing  all  positive  or  all  negative  momentum,  all  but  the  first,  third,  or 
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fourth  terms  may  be  set  to  zero.  Each  of  these  cases  gives  a  simplified  form  of  (83), 
which  can  be  summarized  as 


(r-fi(£))  = 


2k  m 


(  \1/2 

h(IMN) 


vt*  {±k ')  n+  {-^)s±^+k  ] 


(87) 


This  expression  can  now  be  inverted  to  solve  for  the  on-shell  S-matrix  element 

h(MNf  <^-k(£)} 


err  _ 

^ ±k\±k 
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(88) 


For  its  final  fonn,  use  eq.  (78)  to  write  this  as 


err 
° ±k\±k 


h(Av]AV)  "  J  eiEl,hC(t)dt 

_ —co _ 
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(89) 


where  the  correlation  function 

C(t)  =  (i//_\i//+(t))  (90) 

is  defined  between  the  stationary  product  Moller  state  |i//  )  and  the  evolving  reactant 

oo 

Moller  state  |i//,  (0)  •  The  Fourier  transform  of  the  correlation  |  e,E,l^C(t}dt  is 

-oo 

explicitly  a  function  of  energy.  The  remaining  terms  are  implicating  functions  of  energy 
by  the  relation 

k  =  y]2juBii2(E-Er)/U  (91) 

The  internal  energy  E  is  associated  with  any  energy  offset  of  the  surface  as  well  as  any 
internal  degrees  of  freedom. 
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Scattering  on  the  B+H2  Potential  Energy  Surfaces 

The  effective  PES  for  B  +  H2  on  which  the  wavepackets  will  be  propagated  are 

two-dimensional  depending  on  r  and  R.  As  was  used  in  the  derivation  of  the  channel 
packet  method  (see  for  example  eqs.  (70)  and  (71)),  the  product  and  reactant  wavepackets 
are  chosen  to  be  the  direct  product  of  eigenstates  of  the  asymptotic  Hamiltonian.  In  the  r 
dimension  this  is  an  eigenstate  of  the  numerical  potential,  which  must  be  calculated.  As 
expected  from  knowledge  of  the  H2  molecule,  this  eigenstate  will  resemble  the  eigenstate 

of  the  Morse  oscillator.  In  the  R  dimension,  the  eigenstates  are  those  of  the  free  particle. 
As  these  are  not  proper  wavefunctions,  a  linear  combination  is  chosen,  which  for 
mathematical  convenience  will  be  in  the  form  of  a  Gaussian. 

Initially  this  wavepacket  is  chosen  to  be  isolated  to  only  one  PES.  However  eq. 
(43)  shows  that  successive  applications  of  the  full  Hamiltonian  will  mix  the  states,  giving 
nonzero  amplitudes  on  all  PESs.  As  the  wavepacket  is  propagated  forward  in  time,  its 
correlation  function  with  product  states  located  on  each  PES  may  be  calculated 
simultaneously,  allowing  for  the  calculation  of  scattering  matrix  elements  from  the 
reactant  state  to  each  of  the  product  states  with  only  one  propagation  run. 

Propagation  is  facilitated  by  several  practical  considerations,  which  are  addressed 

next. 
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Absorbing  Boundary  Conditions 

Although  the  split-operator  method  of  propagation  takes  advantage  of  the 
relatively  efficient  FFT  algorithm,  it  is  nonetheless  desirable  to  limit  computation  time. 
This  is  accomplished  by  using  the  smallest  practical  grid. 

As  laid  out  in  eq.  (89),  calculating  S-matrix  elements  using  the  CPM  involves 
calculating  the  time-dependent  correlation  function  between  reactant  and  product  Moller 
states.  This  is  typically  done  by  choosing  one  Moller  state  to  evolve  in  time,  while 
leaving  the  other  alone.  Since  the  correlation  function  is  defined  over  all  time,  from  -qo 
to  qo  ,  strictly  speaking  the  evolving  Moller  state  needs  to  be  propagated  over  all  time. 
However,  the  product  Moller  states  may  be  selected  to  be  localized  [24],  so  propagation 
need  only  be  done  over  the  time  when  the  correlation  function  is  nonzero  -  that  is,  until 
the  evolving  Moller  state  has  propagated  sufficiently  far  out  of  the  interaction  region  such 
that  there  is  no  significant  overlap. 

In  many  cases  this  still  requires  propagating  for  a  relatively  long  time.  While 
waiting  for  the  last  residues  of  the  evolving  Moller  state  to  leave  the  interaction  region, 
the  rest  of  the  Moller  state  may  have  diffused  over  a  large  area.  Since  the  split-operator 
propagator  uses  the  FFT,  which  operates  on  a  finite  set  of  points  under  the  assumption 
that  the  complete  data  set  is  periodic  [22],  should  any  portion  of  the  evolving  Moller  state 
reach  the  edge  of  the  grid  it  will  erroneously  reappear  on  the  other  side.  It  will  then 
continue  to  propagate,  potentially  reentering  the  interaction  region  and  introducing  error 
into  the  correlation  function. 

Simply  increasing  the  grid  size  to  accommodate  the  evolving  Moller  state  over  the 
entire  time  interval  necessary  to  calculate  the  correlation  function  will  greatly  increase 
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the  computational  effort  required.  A  more  efficient  solution  is  to  introduce  absorbing 
boundary  conditions  (ABCs)  [24].  ABCs  are  nonphysical  potentials  that  are  added  near 
the  edge  of  the  grid  so  that  when  evolving  states  approach  them  they  are  attenuated. 

ABCs  are  placed  outside  the  interaction  region  in  a  location  where  the  product  Moller 
state  has  zero  amplitude.  They  must  also  be  located  so  that  they  do  no  attenuate  portions 
of  the  wavepacket  which  will  eventually  enter  the  region  of  the  product  Moller  state. 

These  rules  ensure  the  correlation  function  between  the  product  Moller  state  and  evolving 
reactant  Moller  states  is  not  affected,  causing  the  ABC  to  not  alter  the  calculation  of  the 
S -matrix  elements. 

Absorbing  boundary  conditions  are  imaginary-valued,  such  that  the  new  potential 
takes  the  fonn 

V  =  V0±iVabc  (92) 


where  V0  is  the  original  potential.  Adding  the  ABC  causes  the  time-evolution  operator  to 
become 

_  g-iHAt/h 

-if At  /  h-f  ( V0  ±i  Vabc  )  At  /  h 

=  e  v  ' 

The  sign  is  chosen  depending  on  whether  the  propagation  will  go  forward  or  backward  in 
time.  If  At  is  positive,  the  time-evolution  operator  will  attenuate  the  wavepacket  if  Vabc  is 

chosen  to  be  negative.  If  At  is  negative,  Vabc  is  chosen  to  be  positive. 

The  ABC  must  be  chosen  so  that  it  completely  attenuates  the  wavepacket.  Its 
amplitude  must  be  sufficient  to  ensure  it  does  not  transmit  any  of  the  wavepacket,  and  its 
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shape  should  be  such  that  it  does  not  reflect  the  wavepacket.  A  common  choice  is  that  of 
the  Gaussian  -  in  one  dimension  it  would  take  the  form 


Kbc  (*)  =  A? 


-(x-x0)  IB 


(94) 


Where  A  and  B  are  chosen  to  meet  the  criteria  described  above.  If  the  Gaussian  is  too 
narrow,  some  of  the  incoming  wavepacket  will  be  reflected.  If  it  is  too  broad,  the  ABC 
may  intrude  into  the  interaction  region,  introducing  error. 

Other  Practical  Grid  Considerations 

The  use  of  absorbing  boundary  conditions  removes  the  necessity  of  using  a  huge 
grid  by  attenuating  any  part  of  the  wavepacket  that  propagates  well  away  from  the 
interaction  region.  However  there  are  still  several  practical  considerations  when  selecting 
a  grid  for  propagation. 

Mailer  States  and  the  Asymptotic  Limit 

Calculating  S-matrix  elements  using  the  CPM  involves  obtaining  the  time- 
dependent  correlation  function  between  the  reactant  and  product  Moller  states.  Logically 
the  first  step  in  this  process  is  to  calculate  the  Moller  states.  Judicious  selection  of  the 
initial  states  can  greatly  simplify  this  step,  and  will  play  a  small  role  in  grid  selection. 

The  Moller  states  are  defined  (see  eq.  (58)  along  with  the  definition  of  the  Moller 
operators  in  eqs.  (55)  and  (56))  by  propagating  the  reactant  (product)  state  backwards 
(forwards)  in  time  to  t  =  m»  under  the  asymptotic  Hamiltonian  H0 ,  then  forwards 

(backwards)  to  t  =  0  under  the  full  Hamiltonian  H  .  It  is  straightforward  to  perfonn  these 
propagations  using  the  split-operator  propagator.  However  the  required  computation 
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time  (and  the  small  but  theoretically  nonzero  error  associated  with  the  propagator)  may 
be  avoided  entirely  by  defining  the  initial  states  to  already  be  in  the  asymptotic  limit. 

This  choice  means  that  during  the  propagation  H  =  H0  and  the  Moller  operators  are  the 
identity  operator,  and  the  Moller  states  are  simply  the  initial  product  and  reactant  states. 

As  previously  stated,  it  is  desirable  to  limit  the  size  of  the  coordinate  grid  in  order 
to  reduce  the  computation  time  required.  Of  course  it  is  necessary  to  make  the  grid  big 
enough  to  ‘fit  the  problem  on’.  Practically  speaking  this  means  it  must  extend  far  enough 
from  the  interaction  region  to  allow  the  product  and  reactant  states  to  be  defined  in  the 
asymptotic  limit.  The  grid  is  further  extended  beyond  this  point  in  order  to  add  absorbing 
boundary  conditions,  which  must  be  located  where  the  non-evolving  Moller  state  has 
zero  amplitude  in  order  to  not  interfere  with  the  calculation  of  the  correlation  function. 

Fourier  Transform  Pairs:  The  Coordinate  and  Momentum  Grids 

For  each  dimension  the  coordinate  (and  momentum)  grid  are  completely  specified 
by  any  two  of  the  following  parameters  (for  each  dimension  of  the  propagation  space): 
the  range  ( xmax  -xmin),  the  interval  between  adjacent  points  (Ax) ,  and  the  total  number  of 

gridpoints  ( Nx ) .  Given  any  two,  the  third  parameter  is  simply  calculated  by  some  form 
of  the  equation 

x  —  X 

_  max  min 

Nx-1 

The  choice  of  which  parameters  to  specify  is  arbitrary,  although  for  practical 
considerations  directly  specifying  the  range  is  often  advantageous.  Although  it  may  seem 
more  natural  to  specify  the  interval  Ax  and  let  this  detennine  the  number  of  points, 
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frequent  use  of  the  FFT  algorithm  means  it  will  be  advantageous  to  have  Nx  be  a  power 
of  2  [22].  Thus  Nx  is  often  the  second  parameter  specified. 

As  seen  by  eq.  (52),  the  coordinate  and  momentum  grids  are  related  by  the 
Fourier  transfonn.  Therefore  when  selecting  the  coordinate  grid,  the  momentum  grid  is 
defined  as  well.  The  relation  between  the  two  is  given  by  [22] 

A  k  = - — -  (96) 

x  — X 

max  min 

The  number  of  points  in  the  momentum  grid  is  of  course  the  same  as  the  number  of 
points  in  the  coordinate  grid. 

A,  =  Nx  (97) 

These  two  parameters  completely  specify  the  momentum  grid.  The  range  of  momentum 
is  calculated  from  A k  and  Nk  along  with  knowledge  that  the  momentum  grid  will  be 
centered  around  zero: 

(98) 

This  gives  another  piece  to  consider  when  selecting  the  grid.  Since  the  range  of 
momentum  is  never  directly  chosen  and  is  instead  detennined  by  the  fineness  of  the 
coordinate  grid,  care  should  be  taken  to  ensure  that  during  the  propagation  the 
wavepacket  never  reaches  momentum  values  greater  than  that  given  by  eq.  (98).  As 
previously  discussed,  the  FFT  assumes  periodicity,  and  similar  to  the  case  of  the 
wavepacket  going  off  the  grid  in  the  coordinate  representation  which  motivated  use  of 
absorbing  boundary  conditions,  the  wavepacket  can  go  off  the  grid  in  the  momentum 
representation  introducing  error  into  the  propagation.  In  the  case  of  the  coordinate  grid, 
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the  localized  nature  of  the  product  states  meant  portions  of  the  wavepacket  reaching  the 
edge  of  the  grid  were  no  longer  relevant  to  the  calculation  of  the  correlation  function 
between  the  Moller  states,  and  thus  could  be  attenuated  by  employing  absorbing 
boundary  conditions.  No  equivalent  argument  exists  to  allow  us  attenuation  in  the 
momentum  representation  and  the  momentum  grid  simply  must  be  large  enough  to 
contain  the  evolving  wavepacket. 

Fourier  Transform  Pairs:  The  Energy  and  Time  Grids 

Besides  transforming  between  the  coordinate  and  momentum  representations  of 
the  wavepacket  for  use  in  the  split-operator  propagator,  another  use  of  the  Fourier 
transform  is  encountered  when  calculating  S-matrix  elements  using  the  channel  packet 
method.  That  is  the  Fourier  transfonn  of  the  time-dependent  correlation  function  (eq. 

(90)),  which  maps  the  correlation  function,  as  a  function  of  time,  to  the  S-matrix 
elements,  which  are  functions  of  energy.  Thus  time  and  energy  are  Fourier  transform 
pairs,  and  share  the  same  relationship  as  the  coordinate/momentum  pairs: 

A E  =  — — —  (99) 

t  —t 

max  min 

As  previously  discussed,  the  correlation  function  is  defined  over  all  time  but  is 
only  nonzero  when  the  Moller  states  overlap  significantly.  The  time  when  the  correlation 
function  is  nonzero  determines  the  propagation  time  range  (tmax  -  tmin ) .  If  the  initial 

reactant  and  product  states  are  chosen  to  be  located  in  the  asymptotic  limit,  by  definition 
the  correlation  function  is  zero  at  that  time  and  all  of  the  overlap  between  the  fixed  and 
evolving  Moller  state  lies  in  the  future.  In  that  case  tmm  =  0  and  tmax  must  be  sufficiently 
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large  to  capture  all  overlap  between  the  evolving  states.  That  is,  propagation  continues 
until  the  correlation  function  goes  to  zero. 

Care  should  be  taken  to  not  truncate  the  correlation  function  prematurely.  It  is  not 
uncommon  to  encounter  long,  low-amplitude  correlations  as  residuals  of  the  evolving 
wavepacket  oscillate  in  the  interaction  region.  Although  small  in  magnitude  they  recur 
over  a  long  period  of  time,  so  they  contribute  significantly  to  the  integral  over  all  time. 

The  time  range  is  therefore  detennined  by  the  requirement  to  capture  the 
correlation  function  all  times  when  it  is  nonzero.  The  other  piece  needing  to  be  specified 
is  either  the  time  step  A t  or,  equivalently,  the  number  of  points  Nt .  Since  the  accuracy  of 

the  split-operator  depends  on  the  smallness  of  At  (see  eq.  (47)),  the  time  step  is  often  the 
parameter  selected.  The  time  range  and  time  step  completely  specify  the  energy  grid  as 
given  by  eq.  (99)  along  with 

NE=Nt  (100) 

and  knowledge  that  the  energy  grid  will  be  centered  around  zero. 

If  the  time  range  (tmax  -tmm )  is  fixed  by  the  requirement  to  capture  all  of  the 

overlap  between  the  fixed  and  evolving  Moller  state  and  At  is  made  successively  smaller 
in  the  interest  of  greater  accuracy  of  propagation,  the  energy  grid  is  made  successively 
larger.  It  may  be  the  case  that  the  requirement  of  the  split-operator  propagator  for  a  small 
At  will  result  in  the  energy  grid  being  much  larger  than  necessary.  In  that  case,  the 
correlation  function  may  be  downsampled  before  undergoing  the  Fourier  transform  when 
calculating  S-matrix  elements.  Equivalently,  the  correlation  function  need  not  be 
calculated  at  every  propagation  time  step. 
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III.  Results  and  Discussion 


The  B+H2  Effective  Potential  Energy  Surfaces 
Obtaining  the  Surfaces 

The  diabatic  effective  potential  energy  surfaces  Vfff  ( r,R )  as  given  by  eq.  (42)  are 

calculated  by  Fortran  code  developed  by  Garvin  [1],  The  first  part  of  this  work  was  to 
create  a  shell  script  to  interact  with  that  existing  software  to  facilitate  specifying  a 
coordinate  grid  and  basis  size  n  and  to  save  the  resulting  surfaces  in  a  usable  fonnat. 

First  diabatic  surfaces  as  a  function  of  r  and  R  are  calculated  for  every  element  in  the 
n  x  n  diabatic  effective  potential  energy  matrix.  At  every  point  on  the  coordinate  grid 
this  matrix  is  diagonalized.  The  resulting  diagonal  elements  form  the  n  adiabatic 
potential  energy  surfaces  V^f  ( r,R )  which  are  saved  along  with  the  diabatic-to-adiabatic 

transfonnation  U4D  {r,  R)  which  diagonalized  the  potential  matrix  at  that  point. 


The  Diabatic  Surfaces 

/f  ja' 

A  representative  selection  of  diabatic  effective  PESs  (  ,  ( 

shown  here  for  several  values  of  j,k,ja,co.  Both  diagonal  surfaces  and  nondiagonal 
coupling  surfaces  are  considered. 


are 


The  ground  state  surface 


quickly  resolves  to  the  asymptotic 


H2  potential  as  R  increases. 
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Contours  are  shown  for  values  from  0  to  0.27  au,  with  lines  every  0.01  au.  The 
process  of  expanding  the  numerical  surfaces  in  a  truncated  basis  results  in  nonphysical 
features  at  energies  greater  than  0.27  au.  Since  this  energy  is  much  higher  than  those 
accessible  in  the  collision  considered  the  surfaces  have  been  truncated  at  0.27  au.  As  a 
result  in  certain  regions,  for  example  the  region  r  <  3  au ,  R  <  2  au  in  this  plot,  contour 
hatching  is  visible  where  the  surface  has  a  constant  value. 
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The  excited  surface 


1 2  | 
\o’i 


^M) 


2  \\ 
°’i/ 


is  qualitatively  similar  to  the  ground 


state.  It  primarily  differs  in  the  region  of  larger  r,  and  the  asymptotic  II 2  potential  has 


been  raised  in  energy  and  slightly  broadened  due  to  the 


hVCZ  +  l) 

2 II  h/' 


tenn  in  the 


Hamiltonian. 
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Contours  are  shown  for  values  from  0  to  0.27  au,  with  lines  every  0.01  au. 
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Off-diagonal  terms  in  the  diabatic  effective  potential  energy  matrix  represent  coupling 
surfaces.  A  selection  of  them  are  shown  here. 


0.02  — 


Figure  9.  The  Coupling  Surface 
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The  coupling  surface  is  near  zero  for  the  region  R  >  6  au  .  This  indicates  that  a 


wavepacket  initially  located  on  the  ground  state 


will  couple  only  very  weakly  to 


the  state 


in  this  region. 
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Contours  are  shown  from  -0. 12  au  to  .02  au  with  lines  every  0.005  au.  Contour  line 
hatching  is  again  visible  in  regions  where  the  surface  has  a  constant  value. 
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Figure  10.  The  Coupling  Surface 
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Figure  11.  The  Coupling  Surface 
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Contours  are  shown  from  -0.01  au  to  0.025  au  with  lines  every  0.001  au. 
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Figure  12.  The  Coupling  Surface  (  ,  1 

0  "T 


>£M) 


i 

1  -w 


Contour  Plot 
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The  Adiabatic  Surfaces 

While  the  diabatic  surfaces  are  calculated  directly  from  eq.  (42),  the  adiabatic 
surfaces  are  obtained  by  diagonalizing  the  diabatic  effective  potential  energy  matrix  at 
every  point  on  the  coordinate  grid.  This  results  in  several  key  differences  between  the 
diabatic  and  the  adiabatic  surfaces.  Unlike  the  diabatic  surfaces  the  adiabatic  surfaces 
change  as  a  function  of  the  basis  size  n  -  while  increasing  n  simply  results  in  more 
diabatic  surfaces,  the  adiabatic  surfaces  change  as  the  transformation  mixes  in  the 
additional  diabatic  surfaces.  Additionally  there  is  no  clear  labeling  scheme  as  each 

adiabatic  surface  is  a  mixture  of  diabatic  states  with  varying  y  —  ^  basis  states. 

k  co/ 

Often  they  are  ordered  by  increasing  energy  as  this  is  a  common  scheme  for  sorting 
eigenvalues  used  by  diagonalization  algorithms.  In  fact  this  ordering  is  arbitrary,  and  it  is 
sufficient  that  the  adiabatic  potential  energy  matrix  be  diagonal  at  each  point  on  the 
coordinate  grid  and  that  the  unitary  transfonnation  that  diagonalized  the  matrix  at  that 
particular  point  is  known. 
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The  following  series  of  figures  shows  how  the  first  (lowest)  adiabatic  surface 
changes  as  the  basis  size  is  increased.  Contours  are  shown  from  0  to  0.27  au  with  lines 
every  0.01  au. 


Figure  13.  First  Adiabatic  Surface  with  Basis  Size  n  —  2  Contour  Plot 
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Figure  14.  First  Adiabatic  Surface  with  Basis  Size  77  =  8  Contour  Plot 


The  most  interesting  change  here  is  the  height  of  the  surface  in  the  region  around 
f?  =  3  au,  7-  =  2.5  au  .  This  indicates  the  potential  energy  that  must  be  overcome  in  order 
to  access  the  region  of  low  energy  in  the  upper  right  comer.  This  area  represents  the 
reaction  B  +  H2  — »  BH  +  H .  As  the  basis  size  increase  this  barrier  to  reaction  decreases. 
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Figure  15.  First  Adiabatic  Surface  with  Basis  Size  77  =  14  Contour  Plot 


Although  the  surface  does  change  slightly  the  addition  of  states  higher  than  j  =  4 , 
which  corresponds  to  a  basis  size  n  =  14  ,  these  additions  do  not  appreciably  decrease  the 
barrier.  Contour  plots  of  this  surface  with  basis  sizes  n  =  20,  26  appear  nearly  identical. 


61 


2  5  3  3.5  4 

r(au) 


0.05 


Figure  16.  First  Adiabatic  Surface  with  Basis  Size  Yl  =  20  Contour  Plot 
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Figure  17.  First  Adiabatic  Surface  with  Basis  Size  Yl  =  26  Contour  Plot 
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Finally  consider  the  lowest  adiabatic  potential  energy  surfaces  calculated  with  a 


basis  size  n  -  32  .  This  is  the  basis  size  used  in  this  work. 
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Figure  18.  First  Adiabatic  Surface  with  Basis  Size  Tl  —  32  Contour  Plot 


Asymptotic  Eigenstates 

In  preparation  for  propagating  on  these  surfaces  eigenstates  of  the  asymptotic 
Flamiltonian  need  to  be  prepared.  In  the  limit  of  large  R  the  diagonal  diabatic  potential 
energy  surfaces  as  given  by  eq.  (42)  become 
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(101) 


U  Ja 

\k  CO 


VD 

y  eff 


j  jg\_j{j  +  l) 

k  0)  2juHr2 


+  VHi{r)  +  Eel+Eso+Ei 


Off 


The  energy  offset  is  chosen  to  be 


E0fr  Eei  VHl{req) 


placing  the  zero  of  the  energy  between  the  spin-orbit  split  j  -  0  states  at  the  minimum  of 
the  VHj  (r)  potential. 

In  the  asymptotic  limit  the  R  cross-section  of  the  surfaces  follow  the  //2  rotor 
spectrum  with  each  level  being  split  by  spin-orbit  coupling. 


x  10'3 


Figure  19.  Diabatic  Potential  Energy  Surfaces  Cross-section  at  V 


The  asymptotic  potential  in  the  R  dimension  is  that  of  the  free  particle.  When 
selecting  a  wavepacket  for  propagation  any  normalizable  superposition  of  plane  waves  is 
appropriate.  For  convenience  a  Gaussian  is  chosen. 


(102) 
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In  the  asymptotic  limit  the  cross-section  of  the  j  =  0  surface  should  be  the 

VH  (r  )  potential.  When  compared  to  the  VH^  (r)  potential  of  Liu-Siegbahn-Truhlar- 

Horowitz  [25]  [26]  (LSTH),  it  was  noted  that  although  there  was  good  agreement  at  the 
points  included  in  the  original  ab  initio  calculation,  the  interpolation  between  those  points 
differed  to  some  degree.  By  supplying  the  interpolation  routine  with  additional  points 
from  the  LSTH  potential,  this  discrepancy  was  resolved  as  shown  in  the  following  figure. 


Figure  20.  Correction  to  Asymptotic  H 1  Potential 


The  cross-sections  for  values  of  j  higher  than  0  are  offset  by  the  appropriate  rotor 

AJ  +  0 


and  spin  orbit  energies,  and  the  shape  is  slightly  changed  by  the  presence  of  the 


2/6// 


term. 
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Figure  21.  Cross-sections  of  the  Diabatic  PESs  for  j  —  0,2,4  at  Limit  of  Large  R 
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Ligure  22.  Detail  of  Diabatic  PES  Cross-sections  for  j  —  0, 2, 4  Around  r  —  1 .402  au  Showing  Spin-Orbit 

Splitting 
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Unlike  the  case  in  the  ^-direction,  the  eigenstates  in  the  r-direction  must  be 
detennined  numerically.  This  is  done  by  expanding  in  a  basis  of  the  eigenstates  of  the 
simple  hannonic  oscillator. 

In  the  asymptotic  limit  and  at  some  fixed  R  the  Hamiltonian  is  given  by 

*  2 

/V  A.  /V  79 

Ha=T  +  Vr=1f-  +  Vr 
2Rh2 

where  Vr  is  a  the  numerical  cross-section  of  diabatic  potential  energy  surface 

(HUk>=^r+^('-)+£» 

The  kinetic  energy  operator  is  the  same  as  in  the  Hamiltonian  for  the  simple  harmonic 
oscillator,  so  we  can  write  the  Hamiltonian  as 

=k*,-(k.+f) 

Where  Hsho  is  the  Hamiltonian  of  the  simple  hannonic  oscillator 


(103) 


H 


sho 


+—kx2 

2 


The  eigenstates  of  the  simple  harmonic  oscillator  are  the  number  states  |  n)  described  by 
the  eigenvalue  equation  Hsho  \n)  =  Q v(n  + 1  /  2)| nj .  These  eigenstates  have  the 
coordinate  representation 


{x\n)  =  y/n 


(x)  =  (2"nlV^r)  !  (Hn(j3x) 


where  [5  = 


and  'Hn  are  the  Hennite  polynomials. 
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The  treatment  of  the  harmonic  oscillator  eigenstates  is  facilitated  by  defining  the 


raising  and  lowering  operators  a  and  it' 

P  {  ~  .  p 

a  =  —j=  x  +  i  — 
V  2  v 

which  operate  on  the  eigenstates  to  give 


P  ( ~  •  P 

a  =—j=  x-i - 

V2  v  wo 


a\n)  =  yJn\n 


a  77)  =  V/7  +  i  77  + 


The  raising  and  lower  operators  are  related  to  the  coordinate  operator  by 


The  eigenstates  of  the  simple  harmonic  oscillator  are  complete.  Representing  Ha 
in  this  basis  gives 

(77 ' | Ha  1 77.)  -  (77. ’ | Hsho  +{vr-  Vsho ) 1 77) 

=  {n ' I  Hsho  1 77.)  -  (77  ’ |  X2 1 77)  +  (77  ’ |  Vr  1 77) 

=  Qy  77  +  —  e> '  , i —  (77 ' I («+  \n)  +  (n'\V r  \n) 

V  2)  2V2 /T  IV  '  1  1  X  1  1  ' 

The  raising  and  lowering  operators  satisfy  the  relation  it  it 1  -  ait  =  /  where  /  is  the 
identity  operator. 


n'\Ha\n)  =  Uco  77  +  -  <5„,„ 
v  1 7 


2  \[2J3 
+  {n'\Vr\n) 


in'  a' a'  +  a1  a  n^  +  (n'  7?) +  ^77'  aa  7?^ 


Moving  into  the  coordinate  representation  gives 
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I  Ha\n)  = 

Q ) 

n  +  — 

1  a  \  / 

l  2J 

1  k 

2^2(3 


(1  +  2  n) 


S.,. 


1  k 


2^2(3 


(n  + 1)1/2  («  +  2)1/2  Sn,n+2  +  n112  (n  - 1)1/2  Sn,n 


-\y/t,Ax)V,¥n{x)dx 


(104) 


which  gives  a  method  for  numerically  determining  each  matrix  element  (n'\Ha  |  n)  .  This 
result  is  exact  if  an  infinite  number  of  basis  states  are  used.  In  practice  the  basis  must  be 
truncated.  The  result  is  an  m  x  m  matrix  (n'\Ha  |  n)  ,  where  m  is  the  number  of  basis 

states  used.  This  matrix  is  then  diagonalized  and  the  resulting  diagonal  elements  are  the 
eigenvalues  of  the  system  while  the  columns  of  the  unitary  transformation  matrix  used  to 
diagonalize  it  give  the  expansion  coefficients  of  its  eigenvectors  in  the  chosen  basis. 

These  eigenvectors  are  labeled  A"2  (r  ) ,  where  i  indexes  the  surface 

This  work  will  consider  energies  on  the  order  of  the  ground  and  first  excited 
vibrational  state.  The  first  two  eigenvalues  and  their  associated  eigenvectors  were 
calculated  for  the  asymptotic  limit  of  large  R  for  every  potential  energy  surface 
corresponding  to  a  diagonal  element  of  the  effective  potential  energy  matrix.  These 
eigenstates  were  observed  to  converge  when  more  than  10  basis  states  of  the  harmonic 
oscillator  are  used  to  represent  the  Hamiltonian.  To  guarantee  convergence  40  basis 
states  were  used  for  the  calculation. 
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The  ground  and  first  excited  state  are  shown  for  the  cross-section  of  the  ground- 


state  potential  energy  surface 


1  o  { 

\°  1 


V°{r,R  =  30m) 


0  1\ 


Figure  23.  Cross-section  of  First  Diabatic  PES  with  First  Two  Vibrational  Eigenstates 


To  facilitate  propagation  on  a  specific  grid  the  numerically  determined  eigenstates 
(r)  are  sampled  at  the  gridpoints  used  for  the  propagation. 

The  numerically  determined  eigenvalues  are  precisely  the  internal  energies 
referenced  in  equation  (91).  They  incorporate  the  offset  of  each  surface  from  zero  in  the 
asymptotic  limit  as  well  as  the  energy  associated  with  the  internal  vibrational  degree  of 
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freedom.  These  numerically  determined  total  internal  energies  are  shown  in  Table  1  for 


various  values  of  y  - 


j 

k 


and  v . 


Table  1.  Summary  of  Internal  Energies 


V  \ 

°  A 

0  ll 

° 

0  i/ 

2  ^ 

0  1/ 

2 

0  il 

4 

0  il 

4 

0  il 

v  =  0 

0.9884 

0.9957 

1.1499 

1.1573 

1.5212 

1.5287 

v  =  1 

2.8848 

2.8922 

3.0383 

3.0457 

3.3911 

3.3984 

Em  are  in  1 0  2  au  . 

\  r 

v  \ 

6  2\ 

0  ll 

6 

0  V 

8  2\ 

0  ll 

8  A 

0  ll 

10  1\ 

0  ll 

10  f\ 

0  1/ 

v  =  0 

2.0895 

2.0969 

2.8360 

2.8433 

3.7384 

3.7458 

v  =  1 

3.9305 

3.9379 

4.6388 

4.6461 

5.4944 

5.5017 

Since  the  internal  energies  incorporate  the  total  energy  offset  as  well  as  the  energy 
of  the  internal  degrees  of  freedom  they  represent  the  minimum  energy  required  to  access 
those  states.  They  therefore  provide  a  useful  guide  for  determining  ymax  in  order  to 
truncate  the  effective  potential  energy  matrix.  In  this  work  ymax  =  10  is  used,  giving  a 
basis  size  n  =  32  . 
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Propagation 


The  Initial  Wavepacket 


In  preparation  for  propagation  the  initial  wavepacket  is  prepared  in  the  diabatic 
coordinate  representation.  The  selection  of  ymax  =10  means  the  wavepacket  is 
represented  by  a  32  x  1  vector.  This  work  investigates  scattering  matrix  elements  from 


the  state  y  = 


0 

0  i/ 


,  which  will  subsequently  be  referred  to  as  the  ground  state.  The 


initial  wavepacket  is  defined  to  have  amplitude  only  on  this  surface. 


¥-n{r,R)  = 


Y°{r,R)r 

0 

□ 

0 


On  the  ground  state  the  initial  wavepacket  y/°  (r,  R)  is  selected  to  be  an 

eigenstate  of  the  asymptotic  Hamiltonian.  As  previously  discussed,  this  wavepacket  is 
fonned  from  the  direct  product  of  a  Gaussian  in  the  indirection  and  the  ground  state  of 
the  asymptotic  VH  potential  in  the  r-direction: 

An  initial  offset  R0  =  30  au  is  selected  to  place  the  wavepacket  in  the  asymptotic  limit 
along  with  a  spread  parameter  5  =  0.35  and  momentum  offset  k‘(’'  =  -6.75  . 


(105) 


(106) 
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Figure  24.  Initial  Wavepacket  on  the  Ground  State  at  t  —  0  au 


The  Mailer  states 


Since  the  wavepacket  is  located  in  the  asymptotic  limit,  it  is  also  the  reactant 
Moller  state. 


y'-)  =  &-\y'rS)  =  \y'r*) 
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The  product  Moller  states  are  selected  to  have  the  same  form  in  the  /^-dimension  but  with 
a  momentum  offset  kf  =  -k"’  =  +6.75  .  As  with  the  reactant  Moller  states,  their  location 
in  the  asymptotic  limit  eliminates  the  need  to  numerically  propagate  these  states  forwards 
and  backwards  in  time.  Two  product  Moller  states  are  defined  for  every  surface,  one 
with  an  r-dimension  cross-section  of  the  vibrational  ground  state  and  another  with  a 
cross-section  of  the  first  excited  vibrational  state. 

The  width  of  the  momentum  representation  of  the  Moller  states  will  detennine  the 
range  of  energies  considered  by  the  propagation.  Consider  again  the  equation  for 
calculating  scattering  matrix  elements  using  the  channel  packet  method: 


\n  00 

h(|*,.||*r|)  J  e'EllhC(t)dt 

c  fr  — _ _ 

±k’±k  27tmr]-*  (±k')rj+  (±k) 

Here  all  terms  are  shown  explicitly  as  functions  of  energy.  The  denominator  contains  the 
expansion  coefficients  rf  of  the  momentum  representation  of  the  Moller  states  in  the  R- 
dimension.  Since  they  were  chosen  to  be  a  Gaussian  in  the  coordinate  representation, 
these  expansion  coefficients  are  known  -  they  are  simply  given  by  the  inverse  Fourier 
transfonn  of  the  coordinate  Gaussian,  which  is  known  analytically: 


r/±(k)  = 


1/4 


-iR^k 


V  n  J 


To  see  the  range  of  energies  contained  in  these  states  the  if  are  expressed  as  a  function 


of  energy  using  the  relation  k  =  J2jUbh  (E  -  Eml )  /  □. 


(107) 


(108) 
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Figure  25.  PR  Cross-section  of  Initial  Wavepacket 


This  range  of  energies  is  selected  to  be  broad  but  contained  below  the  cutoff  energy 
associated  with  the  truncation  of  the  potential  energy  matrix  at  ymax  =  10. 

Another  consideration  when  selecting  these  paramaters  is  the  requirement  derived 
from  the  channel  packet  method  that  the  wavepacket  momentum  be  either  all  positive  or 
all  negative.  This  requirement  must  be  balanced  against  the  desire  to  obtain  usable  data 
at  low  energies.  The  division  by  rj~*  (±k(fs))  in  eq.  (107)  results  in 

division  by  near-zero  noise  at  the  edge  of  the  Gaussian.  It  is  therefore  desirable  to  locate 
the  Gaussian  such  that  is  has  a  non-negligable  amplitude  at  k  =  0  — >  E  =  Eint . 

The  Reaction  Channel 

While  the  selected  wavepacket  only  has  energies  less  than  those  associated  with 
the  j=12  states,  it  does  contain  energies  higher  than  the  barrier  to  reaction  seen  in  the 
lowest  adiabatic  potential  energy  surface: 
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Figure  26.  Barrier  to  Reaction  in  Lowest  Adiabatic  Potential  Energy  Surface 

The  initial  wavepacket  is  centered  at  r  =  1.402  au.  R  -  30  au  and  is  moving  entirely  in 
the  -R  direction.  As  it  enters  the  interaction  region  portions  of  the  wavepacket  with 
energies  higher  than  □  0.025  au  can  move  into  the  region  of  low  potential  around 
r  >2.5  au,  R  <  2  au  and  exit  the  grid.  Crossing  into  this  region  represents  the  reaction 
B  +  H2  — »  BII  +  H .  Propagation  on  the  BH  +  H  surfaces  is  beyond  the  scope  of  this 
work,  so  the  simplifying  assumption  that  any  wavepacket  exiting  in  that  direction  will 
completely  react  is  required.  The  investigation  of  whether  portions  of  the  wavepacket 
can  cross  onto  the  BH  +  H  surfaces  and  subsequently  return  and  exit  the  grid  along  the 
original  reaction  channel  represents  an  opportunity  for  further  studies.  That  study  would 
require  the  fonnulation  of  the  BH  +  H  effective  potential  energy  surfaces. 
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Absorbing  Boundary  Conditions 


The  assumption  that  the  portions  of  the  wavepacket  that  cross  the  barrier  will 
entirely  exit  is  implimented  through  the  use  of  absorbing  boundary  conditions.  Without 
the  use  of  absorbing  boundary  conditions  any  portion  of  the  wavepacket  leaving  the  grid 
would  erroneously  appear  on  the  opposite  side  due  to  the  FFT’s  assumption  of 
periodicity,  reflect,  and  return  to  the  interaction  region  by  an  unphysical  mechanism.  To 
prevent  this  an  ABC  is  located  in  the  region  /->3au.  The  shape  of  the  ABC  is  a 
Gaussian 

f-w)2 

V,kJ  (r,R)  =  -iA,e  *' 
with  parameters  A1  =  0.025,  Bx  =  1  /  3  . 

With  this  choice  of  parameters  the  ABC  will  only  absorb  portions  of  the 
wavepacket  which  are  reacting.  These  ABCs  are  made  relatively  wide  to  prevent 
reflection  while  completely  absorbing  anything  entering  the  reaction  region. 

Another  ABC  is  located  at  the  edge  of  the  grid  at  large  R.  The  majority  of  the 
wavepacket  will  not  react  and  will  exit  the  interaction  region  along  this  channel.  This 
ABC  prevents  portions  of  the  wavepacket  exiting  in  this  channel  from  crossing  back  to 
R  —  0  ,  reflecting  off  the  potential  there,  and  reentering  the  interaction  region.  This  ABC 
is  a  broader,  shallower  Gaussian  compared  to  Vabc  l  with  the  form 

Kbc.  2(r,R)=-iA2e  "2 

with  the  parameters  A2  =  0.005,  Bx  =  20 . 


(109) 


(110) 
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The  Time  Grid 


An  investigation  of  the  propagation  on  these  surfaces  concluded  that  error 
associated  with  the  split-operator  approximation  converged  to  zero  around  a  time  step  of 
Atc  =  1.5  au  .  For  these  calculations  this  time  step  was  chosen  to  be  At  =  1  au  .  The  total 

propagation  time  of  Tmax  -  4x  106  was  chosen  to  ensure  all  of  the  wavepacket  exited  the 
interaction  region.  This  long  propagation  time  is  necessary  to  capture  the  entire  period 
where  the  correlation  function  is  nonzero. 


78 


A  summary  of  the  propagation  parameters  is  shown  in  Table  2. 


Table  2.  Propagation  Parameters 


Parameter 

Symbol 

Value  (au) 

Coordinate  gridsize 

N* 

256 

K 

32 

Coordinate  grid  range 

r  .  r 

mm  9  max 

0.6702,  4.402 

D  D 

inin 9  inax 

0,  50 

Time  step  size 

At 

1 

Total  propagation  time 

T 

max 

4xl06 

Initial  wavepacket  position 

rO’R0 

1.402,  30 

Initial  wavepacket  momentum 

7  in  lout 

K0 

□6.75 

Initial  wavepacket  spread 

S 

0.35 

Amplitude,  spread  of  ABC  1 

A\,BX 

0.025,  1/3 

Amplitude,  spread  of  ABC  2 

0.005,20 

Masses 

Pb,h1  ’  Mh2 

3109.0,918.6 

Propagation 

Now  all  the  pieces  are  in  place  to  propagate  on  these  surfaces.  Having  defined  the 
initial  wavepacket  on  the  ground  state 
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vZ(r>R>t  =  o) 


0 

□ 

0 


the  correlation  function  between  the  wavepacket  and  each  of  the  64  (32  surfaces  with 
v  =  0, 1 )  product  Moller  states.  This  gives  the  correlation  functions  at  time  t  =  0 . 

To  begin  the  propagation  algorithm  the  diabatic-to-adiabatic  transformtion 
U4D  (r,  R)  is  applied.  At  every  point  on  the  coordiante  grid  u\D  is  a  32  x  32  matrix, 

which  is  multiplied  by  the  32  x  1  diabatic  wavepacket  giving  the  adiabatic  wavepacket  at 
every  point 


yr\r,R) 


V'A(r’R\ 

□ 


This  is  then  multiplied  by  the  first  of  the  potential  energy  operators,  e  ,v‘ifAn2  which  is 
diagonal  in  this  representation.  This  intennediate  wavepacket  is  then  transfonned  back 
into  the  diabatic  representation.  Next  it  is  moved  into  the  momentum  representation  by 
taking  the  two-dimensional  Fourier  transfonn  of  each  of  its  32  rows  y/A  (r,  R  )  .  In  the 

diabatic  momentum  representation  the  kinetic  energy  operator  is  diagonal. 

e-™/a(kr,kR)  = 

After  transforming  back  to  the  adiabatic  coordinate  representation  the  second  potential 
energy  operator  is  applied  in  the  same  way,  giving  the  evolved  wavepacket  at  the  next 


_i  {kr  +kg/2/jB  Hl^At 


0 


□ 


0 


i^kr  12  +kR  12  Mb,H2 
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time  step  in  the  adiabatic  representation.  This  is  transfonned  to  the  diabatic 
representation  giving  y/D  ( r,R,t  =  At) .  The  correlation  functions  between  this  evolved 

wavepacket  and  each  of  the  product  Moller  states  is  calculated  for  t  -  At .  This  process  is 
iterated  for  the  desired  number  of  time  steps. 

Visualizing  the  Propagation 


The  wavepacket  is  defined  to  be  located  only  on  the  ground  state  at  t  =  0 .  The 
off-diagonal  elements  of  the  diabatic  effective  potential  energy  matrix  have  the  effect  of 
mixing  the  wavepacket  onto  the  other  surfaces.  This  process  is  shown  considering  the 


/o  i  ,  0  -A 

wavepacket’ s  evolution  on  ground  state  diabatic  surface  (  {  Kff  y,R)  l  )  and  the 

W  2  ”  2/ 


/o  I  ,  0  f\ 

next  higher  surface  (  :  Vej]  R)  {  )  associated  with  the  boron  electronic  fine 

W  2  0  2 


structure  transition.  At  every  time  point  shown  the  correlation  function  between  the 


evolving  reactant  state  and  the  stationary  product  Moller  state  (located  at  R  —  30  au  with 


momentum  in  the  positive  R  direction)  is  shown  from  time  t  =  0  au  up  to  the  time  being 
considered.  The  product  Moller  states  are  represented  by  a  single  black  contour. 
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At  t  =  0  the  wavepacket  is  placed  in  the  asymptotic  limit  and  is  located  only  on 
the  ground  state.  The  initial  correlation  is  zero  because  the  wavepacket  has  only  negative 
momentum  while  the  product  Moller  state  has  only  positive  momentum. 


Figure  27.  Wavepacket  on  the  First  Two  Surfaces  at  t  —  0  au 


82 


At  t  =  5000  au  the  wavepacket  has  moved  towards  the  interaction  region.  No 
significant  coupling  to  higher  surfaces  or  correlation  with  the  product  Moller  states  is 
seen. 


Figure  28.  Wavepacket  on  the  First  Two  Surfaces  at  t  =  5000  au 
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At  t  =  12500  au  the  wavepacket  has  entered  the  interaction  region  and  has  coupled 


to  the  higher  surfaces. 


3  4 

t(au) 


1  2  3  4  5  6  7 

I  (au)  x  10« 


Figure  29.  Wavepacket  on  the  First  Two  Surfaces  at  t  —  12500  au 


84 


For  t  =  20000  -  45000  au  the  wavepacket  is  exiting  the  interaction  region  on  all 
of  the  surfaces.  A  significant  portion  is  located  in  the  region  of  the  product  Moller  states 
resulting  in  the  most  significant  amplitude  in  the  correlation  functions. 


Figure  30.  Wavepacket  on  the  First  Two  Surfaces  at  t  =  20000  au 
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Figure  31.  Wavepacket  on  the  First  Two  Surfaces  at  t  =  27500  au 
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Figure  32.  Wavepacket  on  the  First  Two  Surfaces  at  t  —  35000  au 
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Figure  33.  Wavepacket  on  the  First  Two  Surfaces  at  t  =  42500  au 
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At  t  =  50000  au  most  of  the  wavepacket  has  passed  by  the  product  Moller  states 


and  the  correlation  function  returns  to  low  amplitude. 


Figure  34.  Wavepacket  on  the  First  Two  Surfaces  at  t  —  50000  au 


89 


At  t  -  65000  au  the  wavepacket  is  no  longer  visible  on  this  scale.  However  the 
remaining  low-amplitude  long-duration  features  of  correlation  function  play  an  important 
part  in  forming  the  scattering  matrix  elements. 


Figure  35.  Wavepacket  on  the  First  Two  Surfaces  at  t  =  65000  au 
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The  Scattering  Matrix  Elements 


Calculating  the  Scattering  Matrix  Elements 

The  correlation  functions  between  the  product  Moller  states  and  the  evolving 
reactant  Moller  states  is  used  to  calculate  the  scattering  matrix  elements  using  the  relation 


•s214(£)  =  - 


h(|*r,||*7|)'':  |  e'E""C(t)dt 


Innir]  * (±k'^rj+  (±k) 

First  the  Fourier  transform  of  each  correlation  function  is  taken.  These  are  multiplied  by 
the  appropriate  prefactor  according  to  the  expansion  coefficients  of  the  respective 


reactant  and  product  Moller  states,  rf  [kR )  = 


/  .  \  l/4 

f  28 2A 


S2  f  j  iRVlkR 


.  These  are 


V  n  J 


implicitly  functions  of  energy  using  the  relation  k  =  (  E  -  Eml )  /  □  where  the 


internal  energies  were  calculated  numerically  and  are  summarized  in  Table  1. 


The  Scattering  Matrix  Elements  for  B+Hj 


The  scattering  matrix  elements  represent  probability  amplitudes  for  scattering  to  a 


given  state.  Here  they  represent  transitions  from  the  ground  state 


J  Ja 

k  o) 


,v, 


0  1 
0  i 


,0 


to  each  of  the 


J  Ja 

k  03 


,v)  states.  The  probability  of  these  transitions  is  shown  by  the 


amplitude  squared  of  the  scattering  matrix  element  and  are  obtained  as  a  function  of 


energy. 
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Figure  40.  Probability  to  Transition  from 
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Figure  41.  Probability  to  Transition  from 


0 

0 


10 


to  the  Six  Allowed 


j 


States 


These  scattering  matrix  elements  are  qualitatively  similar  to  previous  calculations 
using  a  one  dimensional  propagator  with  the  restriction  r  =  req  ,  with  broad  Stueckelberg- 
type  [27]  oscillatory  behavior  as  a  function  of  energy,  along  with  high-frequency 
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Feshbach  resonances  [28]  corresponding  to  portions  of  the  wavepacket  which  couple  to 
higher  energy  surface  and  become  trapped  in  quasi-stable  states.  These  Feshbach 
resonances  are  the  primary  driver  for  the  long  total  propagation  time,  and  are  discussed  in 
greater  detail  in  a  subsequent  section. 

Convergence  of  Basis  Size 

The  effective  diabatic  potential  energy  matrix  is  infinite  in  dimension.  Energetic 
considerations  motivated  the  truncation  of  this  matrix  at  some  /  .For  this  work  /  is 

J  max  J  max 

chosen  to  be  10,  giving  a  basis  size  n  =  2  +  3 ymax  =  32  .  Repeating  the  calculation  of  the 

scattering  matrix  elements  for  bases  of  increasing  size  demonstrates  that  this  truncation  is 
appropriate  over  specific  energy  ranges. 

Consider  the  scattering  matrix  element  for  a  wavepacket  entering  and  leaving  on 
the  ground  state.  This  case  is  often  referred  to  as  reflection.  The  probability  for 
reflection  as  a  function  of  energy  is  shown  for  calculations  using  /max  =  0, 2, 4, 6, 8, 10 , 
which  result  in  basis  sizes  n  =  2,8,14,20,26,32  respectively. 
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Figure  42.  Comparison  of  Reflection  Probability  for  =0,2,4,6,8,10 


max 


The  energy  0.009883  au  represents  the  internal  energy  of  the  ground  state,  which 
will  serve  as  a  constant  offset  from  zero.  The  first  line  to  diverge  is  that  of  the 
calculation  using  ymax  =  0  .  As  expected  this  occurs  near  the  energy  where  the  j  =  2  states 

become  accessible.  As  energy  increases  calculations  based  on  lower  numbers  of  basis 
states  diverge  indicating  the  energy  range  over  which  they  are  valid  has  been  exceeded. 
At  no  point  in  this  range  do  the  ymax  =  8  values  diverge  from  the  /max  =  10  values, 

indicating  that  the  calculation  has  converged  over  this  energy  range. 

Examination  of  the  scattering  matrix  elements  for  the  fine-structure  transition  to 
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Figure  43.  Comparison  of  Fine  Structure  Probability  for  /max  =  0, 2, 4, 6,8,10 


Sum  of  Reaction  Probabilities 

To  the  extent  that  all  possible  reactions  have  been  accounted  for  the  sum  of  the 
amplitude  squared  of  all  scattering  matrix  elements  should  be  1 .  This  is  a  useful  check 
for  consistency  and  may  be  used  as  a  diagnostic  tool  to  check  for  convergence  based  on 
time  step  size  and  total  propagation  time.  Prematurely  stopping  the  propagation  truncates 
the  correlation  functions,  and  their  Fourier  transform,  and  subsquently  the  scattering 
matrix  elements,  will  exibit  ringing.  Selection  of  an  inappropriatly  large  time  step 
introduces  an  erroneous  phase  shift  in  the  correlation  funtion.  This  results  in  a  translation 
occuring  to  the  Fourier  transform.  While  this  shift  may  not  be  noticible  when  examing 
individual  scattering  matrix  elements,  the  requirement  that  they  sum  to  unity  makes  it 
apparent.  The  equation  for  the  scattering  matrix  elements 

oo 

0|k(fs)|  |  e'Et  (i//_  \y/+{t)) dt 
S±k'±k  ^  2^^  (±k\E))ri+  (±*(£)) 
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involves  dividing  the  curve  associated  with  the  Fourier  transform  of  the  correlation 
function  by  the  product  of  two  Gaussians,  so  it  is  sensitive  to  translation  of  the  numerator 
with  respect  to  the  denominator  -  instead  of  summing  to  one  the  scattering  matrix 
elements  will  curve  through  it. 

The  sum  of  the  scattering  matrix  elements  for  this  calculation  is  shown  in  Figure 

44. 


Figure  44.  Sum  of  Scattering  Matrix  Elements 

High  frequency  oscillations  are  visible  at  the  internal  energy  of  the  ground  state 
E  associated  with  division  by  near-zero  error  where  the  Gaussian  in  the  denominator  of 
eq.  XX  has  reached  a  very  small  value.  This  situation  is  repeated  at  energies  for  the 
higher  states  E j=2  and  F;  4 .  Otherwise  the  sum  is  equal  to  one  up  to  a  threshold  energy 
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associated  with  the  wavepacket  overcoming  the  barrier  to  reaction  and  being  attenuated 
by  the  absorbing  boundary  condition.  Since  this  work  does  not  consider  the  BH  +  H 
surfaces,  it  is  not  possible  to  create  a  product  Moller  state  associated  with  reaction,  which 
necessitates  use  of  ABCs  resulting  in  loss  of  probability  represented  in  the  final  sum. 
Under  the  assumption  that  all  of  the  wavepacket  that  enters  the  reaction  well  will  react, 
the  deviation  of  this  sum  from  one  represents  the  probability  for  reaction. 


Figure  45.  Probability  to  React  to  Form  BH+H 


This  analysis  allows  for  the  prediction  that  reaction  to  BH  +  H  will  begin  to  occur 
with  a  minimum  energy  of  0.022  au. 

Vibrational  Transitions 

As  expected  transitions  to  higher  vibrational  eigenstates  are  not  available  for 
energies  less  than  the  total  internal  energies  associated  with  those  states.  Above 
E  =  0.03  au  these  states  become  available. 
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Figure  46.  Probability  to  Transition  from 
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Figure  47.  Probability  to  Transition  from 


0 

ift\ 

,  ,  0  )  to  the  Six  Allowed 

2 

Ja,l) 

0 

i  / 

j 

CO  j 

100 


Figure  48.  Probability  to  Transition  from 
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Figure  49.  Probability  to  Transition  from 
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As  expected  based  on  their  internal  energies,  there  are  no  transitions  to  the  j  =  8,  v  =  1  or 
j  =  10,  v  =  1  in  this  energy  range.  The  energies  associated  with  transitions  to  the  higher 
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vibrational  eigenstates  are  at  the  upper  limit  of  the  energies  present  in  the  initial 
wavepacket,  so  usable  data  is  obtained  over  a  smaller  range  of  energies. 

Feshbach  Resonances 

Portions  of  the  wavepacket  that  are  coupled  to  higher  surfaces  have  energies  near 
the  eigenvalues  of  the  shallow  wells  on  those  surfaces.  Figure  50  shows  the  wells  for  the 
lowest  three  adiabatic  potential  energy  surfaces  at  r  =  1 .402  au  . 


Figure  50.  Cross-section  of  First  Three  Adiabatic  Surfaces  at  r  —  1 .402  au  Showing  Shallow  Wells 


These  shallow  well  potentials  cause  portions  of  the  wavepacket  to  be  trapped  in  a 
quasi-stable  state  and  remain  there  until  they  are  coupled  back  down  to  the  lower  energy 
surfaces  at  which  point  they  will  exit  the  interaction  region.  This  phenomenon  results  in 
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high  amplitude  peaks  appearing  in  the  scattering  matrix  elements  for  those  lower  surfaces 
at  energies  near  the  eigenstates  of  the  wells  of  the  higher  surfaces. 


Energy  (au) 


Figure  51.  Detail  of  Feshbach  Resonances  in  Reflection  Scattering  Matrix  Element 


These  resonances  span  the  energies  of  the 
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and  are  offset  from  the  internal  energies  of  those  states  by  a  small  amount  associated  with 
the  energy  of  the  shallow  potential  wells. 


Although  smaller  in  magnitude,  the 
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Figure  52.  Feshbach  Resonance  Detail  for  Transition 
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Feshbach  resonances  are  also  observed  in  calculations  utilizing  a  one-dimensional 
propagator.  By  including  the  r  dimension,  the  resonances  now  occur  at  energies  that  are 
the  eigenvalues  of  a  two-dimensional  well.  In  the  region  of  the  shallow  potential  wells  in 
the  f? -dimension  the  r  cross-section  is  approximately  that  of  the  asymptotic  1 1 2  potential. 
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Figure  53.  Cross-section  of  First  Three  Adiabatic  Surfaces  at  R  =  7  au 

The  steepness  of  the  r  cross-section  (this  may  be  quantified,  for  example,  by 
representing  it  as  a  simple  harmonic  oscillator  and  examining  the  spring  constant)  is  a 
factor  of  :  104  higher  than  in  the  ^-direction.  Because  of  this  the  two-dimensional  well 
does  not  offer  more  quasi-stable  states  -  they  are  simply  offset  from  zero  by 
approximately  the  energy  of  the  vibrational  ground  state.  Because  of  this  relative 
steepness,  and  due  to  the  fact  that  in  the  region  of  the  shallow  wells  in  the  /^-direction,  the 
shape  of  the  Feshbach  resonances  is  not  significantly  different  in  the  two-dimensional 
case  compared  to  the  one-dimensional  case  obtained  by  constraining  the  molecular 
bondlength  to  r  =  1 .402  au  . 

Because  they  have  distinctive  features  and  are  associated  with  portions  of  the 
wavepacket  that  become  trapped  and  exit  the  interaction  region  slowly,  the  Feshbach 
resonances  provide  a  useful  tool  for  determining  convergence  based  on  total  propagation 
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time.  Figure  54  shows  the  Feshbach  resonances  observed  in  the  scattering  matrix 
element  for  reflection  on  the  ground  state  as  calculated  using  various  total  propatation 


times  T 


Figure  54.  Detail  of  Feshbach  Resonances  for  Tmm  —2,3,4  xlO6  au 


Around  the  sharp  peaks  of  the  Feshbach  resonances,  the  line  corresponding  to 
Tmx  =2x1  ()'’  au  deviates  from  the  other  lines.  Because  the  correlation  function  has  a 
shorter  total  time  its  Fourier  transfonn  does  not  have  the  resolution  necessary  to  capture 
these  features.  The  agreement  between  the  calculation  using  Tmx  =3x1  O'’  au  and 

Tmx  =  4x  106  au  indicates  the  calculation  has  converged  at  that  total  propagation  time. 

Comparison  to  One  Dimensional  Calculation 

Previous  calculations  of  the  B  +  H2  scattering  matrix  elements  were  done  by 
fixing  the  molecular  bondlength  to  its  equilibrium  value  r  —  1 .402  au  and  propagating  on 
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the  resulting  one  dimensional  potential  energy  surfaces.  Including  the  r  degree  of 
freedom  changes  the  scattering  matrix  elements  in  several  ways.  Figure  55  shows  a 
comparison  between  the  one  dimensional  and  two  dimensional  results  for  the  transitions 


from 


0  }  \ 

,0  to 

0  i  / 


0  \ 


0 


,  0  )  and 


0  f 
0  i 


,01.  In  these  plots  the  two  dimensional  results  are 


shifted  by  the  energy  of  the  ground  vibrational  state  in  order  to  make  the  energy  scales 
the  same. 
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Figure  56.  1-D  and  2-D  Results  for  the  Transition  from 


0 

0 


0 

0 


The  scattering  matrix  elements  for  the  one  dimensional  and  two  dimensional 
calculations  are  very  similar  with  the  primary  difference  being  a  shift  towards  lower 
energies  in  the  two  dimensional  case.  This  shift  grows  larger  at  higher  energies. 

h2y'(y  +  i) 

The  two  dimensional  potential  energy  surfaces  contain  a - - — term,  as  seen 

2/A/ 

in  eq.  (3 1).  Besides  providing  the  rotor  energy  separation  to  the  surfaces,  this  term  also 
flattens  out  the  r-cross  section  of  the  higher  surfaces.  As  a  result  of  this  flattening,  the 
ground  vibrational  energy  is  slightly  lowered  on  surfaces  with  higher  j.  For  example,  the 
ground  vibrational  energy  on  the  j  =  2  surfaces  is  :  5  x  1 0  5  au  lower  than  the  ground 
vibrational  energy  on  the  j  =  0  surfaces.  This  effectively  makes  the  upper  surfaces 
accessible  at  lower  energies,  resulting  in  a  shift  in  the  scattering  matrix  elements.  This 
shift  is  examined  at  the  first  set  of  Feshbach  resonances  in  the  ground  state  transition. 
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The  Feshbach  resonances  from  the  two  dimensional  calculation  are  shifted  in 
energy  by  :  5  x  1 0  au  as  predicted.  While  their  structure  is  qualitatively  the  same,  the 
two  dimensional  calculation  yields  additional  low  amplitude  features.  This  is  also  seen  in 
the  Feshbach  resonances  in  the  fine  structure  transition. 


109 


110 


IV.  Conclusion 


Coupling  between  different  degrees  of  freedom  presents  a  challenge  when 
modeling  the  dynamics  of  atomic  and  molecular  collisions.  The  collision  of  boron  and 
molecular  hydrogen  is  an  example  of  a  system  in  which  this  nonadiabatic  behavior  plays 
a  key  role.  This  work  models  the  dynamics  of  B  +  H2  as  they  collide  and  calculates  for  a 

range  of  energies  the  probability  of  changes  to  the  molecule’s  rotational  and  vibrational 
state  as  well  as  the  electronic  fine-structure  transition  of  the  single  open  shell  electron  in 
boron. 

Garvin  [1]  obtained  the  effective  potential  energy  surfaces  for  B  +  IB  by 
combining  the  system’s  rotational  energy  with  the  adiabatic  electronic  potential  energy 
surfaces  calculated  by  Yarkony  and  representing  the  result  in  an  angular  momentum 
basis.  These  effective  PES  are  two-dimensional  in  r,  the  //2bondlcngth,  and  R,  the 

distance  from  the  boron  atom  to  the  H2  center  of  mass.  As  a  result  of  the  angular 

momentum  basis  chosen  these  surfaces  are  labeled  by  the  angular  momentum  quantum 
numbers.  The  PES  are  coupled,  allowing  a  wavepacket  initially  on  a  single  surfaces  to 
transition  to  multiple  surfaces.  Since  there  are  an  infinite  number  of  rotational  states 
there  are  an  infinite  number  of  coupled  PES.  For  practical  considerations  only  those 
surfaces  which  are  energetically  accessible  are  considered.  In  this  work  only  values  of 
total  angular  momentum  J  -  1  /  2  are  considered,  and  the  projection  of  J  onto  the  body 
fixed  z-axis,  labeled  P,  is  assumed  to  be  constant  under  the  centrifugal  sudden 
approximation. 
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Wavepackets  are  propagated  on  these  surfaces  using  a  split-operator  propagator 
which  approximates  the  time-evolution  operator  for  small  time  steps.  Initially  a 
wavepacket  is  located  on  a  single  PES.  This  wavepacket  is  chosen  to  be  an  eigenstate  of 
the  asymptotic  Hamiltonian  and  is  referred  to  as  the  reactant  Moller  state.  As  it  is 
propagated  forward  in  time  it  enters  the  interaction  region  and  is  coupled  to  the  other 
surfaces  so  that  it  then  exits  the  interaction  region  on  all  of  the  surfaces.  At  every  point 
in  time  the  correlation  function  is  calculated  between  the  evolving  state  and  a  product 
state  located  in  the  asymptotic  limit  of  each  PES.  The  Fourier  transfonn  of  these 
correlation  functions  can  be  used  to  calculate  scattering  matrix  elements  using  the 
channel  packet  method.  The  absolute  value  squared  of  the  scattering  matrix  elements  is 
the  probability  that  as  a  result  of  the  collision  the  system  will  transition  from  the  state 
described  by  the  reactant  wavepacket  to  a  state  described  by  the  product  wavepacket. 


Table  3.  Summary  of  Steps  for  Calculating  Scattering  Matrix  Elements 


1. 

Obtain  electronic  adiabatic  PES  for  B  +  1 1 2  (Yarkony) 

2. 

Obtain  effective  diabatic  PES  for  B  +  //,  (Garvin) 

3. 

Numerically  determine  eigenstates  and  eigenvalues  for  each  PES  in  asymptotic  limit 

4. 

Define  reactant  Moller  state  on  ground  state  and  product  Moller  states  (representing  both  the 

v  =  0,1  vibrational  levels)  on  each  PES 

5. 

Propagate  reactant  Moller  state  forward  in  time 

6. 

Calculate  the  correlation  function  between  the  reactant  Moller  state  and  each  of  the  product 

Moller  states  (over  all  time  for  which  it  is  nonzero) 

7. 

Calculate  the  scattering  matrix  elements  using  the  channel  packet  method  (eq.  (89)) 
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The  scattering  matrix  elements  exhibit  broad  oscillatory  behavior  with  rapidly 
oscillating  Feshbach  resonances  associated  with  energies  of  the  quasi-stable  states  of  the 
upper  surfaces.  Compared  to  a  calculation  requiring  r  —  1 .402  au  the  scattering  matrix 
elements  are  shifted  towards  lower  energies  consistent  with  the  lowering  of  the 
vibrational  energies  as  the  rotational  energy  increases.  Including  the  r  degree  of  freedom 
also  has  the  important  result  of  allowing  the  reaction  to  BH  +  H  ,  which  is  observed  for 
energies  beginning  with  0.022  au.  Finally,  considering  the  r  degree  of  freedom  allows 
for  calculation  of  scattering  matrix  elements  for  transitions  involving  changes  to  the 
vibrational  eigenstate. 


Summary  of  Key  Contributions 

This  work  extends  the  channel  packet  method  to  include  treatment  of 
nonadiabatic  systems  in  which  the  reaction  dynamics  are  described  by  potential  energy 
surfaces  that  are  two  dimensional.  This  results  in  a  calculation  based  on  fewer 
unphysical  constraints  and  allows  for  the  consideration  of  changes  of  eigenstates 
associated  with  new  degrees  of  freedom,  as  well  as  capturing  important  reaction 
dynamics  which  may  not  otherwise  be  predicted. 

This  work  implements  this  methodology  for  the  collision  of  B  +  H2,  resulting  in 
transition  probabilities  for  changes  to  the  boron  electronic  fine  structure  and  the  rotational 
and  vibrational  eigenstates  of  the  H2  molecule  as  a  result  of  the  collision.  It  also  allows 
a  prediction  that  reaction  B  +  H^  — »  BH  +  H  will  begin  to  occur  at  a  total  energy  of 
0.022  au . 
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Recommendations  for  Future  Work 


By  extending  the  effective  potential  energy  surfaces  to  include  values  of  total 
angular  momentum  J  greater  than  1/2  total  scattering  cross-sections  may  be  calculated 
and  used  to  predict  reaction  rates. 

The  potential  energy  surfaces  could  also  be  extended  to  include  states  for  which 
there  is  a  change  in  the  projection  of  the  total  angular  momentum  onto  the  body  fixed  z- 
axis  ( P ).  This  would  allow  the  computation  to  be  done  without  invoking  the  centrifugal 
sudden  approximation. 

Calculating  the  BH  +  H  potential  energy  surfaces  would  allow  propagation  on 
these  surfaces  and  allow  for  calculation  of  scattering  events  at  higher  energies  without 
requiring  the  assumption  that  any  portion  of  the  wavepacket  that  crosses  the  barrier  to 
reaction  seen  on  the  lowest  adiabatic  PES  will  react  completely. 
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